Как я могу использовать fsolve для построения решений функции?

Во-первых, вам нужно основное уравнение как функция веса и t, для которого вы хотите fsolve найти ноль. Затем для каждого веса вы можете зафиксировать его в другой функции, которую вы затем решаете для t:

alpha = 0.7e-7;
rbya = 0.464;
k = 1:22;
a = @(weight)(weight./(1360*pi)).^(1/3);
eqn = @(weight,t)2/pi/rbya*sum((-1).^(k-1)./k.*sin(k*pi*rbya).*exp(-1.*k.^2.*pi^2.*alpha.*t./(a(weight).^2)))-57/80;

weights = 4:8;
ts = zeros(size(weights));
for i = 1:numel(weights)
    sub_eqn = @(t)eqn(weights(i),t);
    ts(i)=fsolve(sub_eqn,14400);
end

plot(weights,ts/(60*60))
xlabel("Weight (kg)")
ylabel("Cooking Time (hrs)")

Если вы хотите решить сразу весь набор уравнений, то вам нужно быть осторожным с размерами массивов (как вы уже поняли, подробнее здесь). k должен быть вектором-столбцом, чтобы sum суммировался по каждому столбцу, а weights должен быть вектором-строкой, чтобы поэлементные операции повторяли k для каждого веса. Вам также нужно, чтобы ваш список начальных предположений был того же размера, что и weights, чтобы fsolve мог иметь предположение для каждого веса:

alpha = 0.7e-7;
rbya = 0.464;
k = (1:22)';

a = @(weight)(weight./(1360*pi)).^(1/3);

weights = 4:8;
eqn = @(t)2/pi/rbya*sum((-1).^(k-1)./k.*sin(k*pi*rbya).*exp(-1.*k.^2.*pi^2.*alpha.*t./(a(weights).^2)))-57/80;

ts=fsolve(eqn,repmat(14400,size(weights)));

plot(weights,ts/(60*60))
xlabel("Weight (kg)")
ylabel("Cooking Time (hrs)")

Обратите внимание, что при использовании двух методов вы получаете несколько разные ответы.