How can I plot this loop?

I'm having trouble plotting this loop in Matlab. By doing the fprintf I'm able to make a file with the data but I would like to have Matlab graph it as well. Thanks.

for j=1:nt

      ha=hac(h,rhoa,g,ta,ti,rmua,ka,cpa);
      tf=tfire(t,to);

      k1s=dt*rhss(t,ts,ti,ta,ks,rhos,cps,ls,li,tf);
      k1i=dt*rhsi(t,ts,ti,ta,kbar,rhoi,cpi,ls,li,ha);
      k1a= dt*rhsa(t,ti,ta,rhoa,cpa,la,ha);

      ha=hac(h,rhoa,g,ta+0.5*k1a,ti+0.5*k1i,rmua,ka,cpa);
      tf=tfire(t+0.5*dt,to);

      k2s=dt*rhss(t+0.5*dt,ts+0.5*k1s,ti+0.5*k1i,ta+0.5*k1a,ks,rhos,cps,ls,li,tf);
      k2i=dt*rhsi(t+0.5*dt,ts+0.5*k1s,ti+0.5*k1i,ta+0.5*k1a,kbar,rhoi,cpi,ls,li,ha);
      k2a=dt*rhsa(t+0.5*dt,ti+0.5*k1i,ta+0.5*k1a,rhoa,cpa,la,ha);

      k3s=dt*rhss(t+0.5*dt,ts+0.5*k2s,ti+0.5*k2i,ta+0.5*k2a,ks,rhos,cps,ls,li,tf);
      k3i=dt*rhsi(t+0.5*dt,ts+0.5*k2s,ti+0.5*k2i,ta+0.5*k2a,kbar,rhoi,cpi,ls,li,ha);
      k3a=dt*rhsa(t+0.5*dt,ti+0.5*k2i,ta+0.5*k2a,rhoa,cpa,la,ha);

      ha=hac(h,rhoa,g,ta+k3a,ti+k3i,rmua,ka,cpa);
      tf=tfire(t+dt,to);

      k4s=dt*rhss(t+dt,ts+k3s,ti+k3i,ta+k3a,ks,rhos,cps,ls,li,tf);
      k4i=dt*rhsi(t+dt,ts+k3s,ti+k3i,ta+k3a,kbar,rhoi,cpi,ls,li,ha);
      k4a=dt*rhsa(t+dt,ti+k3i,ta+k3a,rhoa,cpa,la,ha);

      ts=ts+1.0/6.0*k1s+1.0d0/3.0d0*k2s+1.0/3.0*k3s+1.0/6.0*k4s;
      ti=ti+1.0/6.0*k1i+1.0d0/3.0d0*k2i+1.0/3.0*k3i+1.0/6.0*k4i;
      ta=ta+1.0/6.0*k1a+1.0d0/3.0d0*k2a+1.0/3.0*k3a+1.0/6.0*k4a;

      t=t+dt;

      fprintf(fid,'%10.6f %10.6f %10.6f %10.6f %10.6f\r\n',t,tf,ts,ti,ta);

end

[Rik]

If you store the results in vector elements, plotting is trivial. What variable would you like to plot?

[Carl Hendrickson]

I would like to plot t(time) vs t(temp)f, ts, ti, and ta.

Answer 1

[Rik]

There are better ways to do this, but without re-writing a sizable chunk of your code that's a bit difficult. So in the example below I solved it the ugly way. I can't run it as an example, as I'm missing several parameters.

outputs=zeros(nt,5);
for j=1:nt
      ha=hac(h,rhoa,g,ta,ti,rmua,ka,cpa);
      tf=tfire(t,to);
      k1s=dt*rhss(t,ts,ti,ta,ks,rhos,cps,ls,li,tf);
      k1i=dt*rhsi(t,ts,ti,ta,kbar,rhoi,cpi,ls,li,ha);
      k1a= dt*rhsa(t,ti,ta,rhoa,cpa,la,ha);
      ha=hac(h,rhoa,g,ta+0.5*k1a,ti+0.5*k1i,rmua,ka,cpa);
      tf=tfire(t+0.5*dt,to);
      k2s=dt*rhss(t+0.5*dt,ts+0.5*k1s,ti+0.5*k1i,ta+0.5*k1a,ks,rhos,cps,ls,li,tf);
      k2i=dt*rhsi(t+0.5*dt,ts+0.5*k1s,ti+0.5*k1i,ta+0.5*k1a,kbar,rhoi,cpi,ls,li,ha);
      k2a=dt*rhsa(t+0.5*dt,ti+0.5*k1i,ta+0.5*k1a,rhoa,cpa,la,ha);
      k3s=dt*rhss(t+0.5*dt,ts+0.5*k2s,ti+0.5*k2i,ta+0.5*k2a,ks,rhos,cps,ls,li,tf);
      k3i=dt*rhsi(t+0.5*dt,ts+0.5*k2s,ti+0.5*k2i,ta+0.5*k2a,kbar,rhoi,cpi,ls,li,ha);
      k3a=dt*rhsa(t+0.5*dt,ti+0.5*k2i,ta+0.5*k2a,rhoa,cpa,la,ha);
      ha=hac(h,rhoa,g,ta+k3a,ti+k3i,rmua,ka,cpa);
      tf=tfire(t+dt,to);
      k4s=dt*rhss(t+dt,ts+k3s,ti+k3i,ta+k3a,ks,rhos,cps,ls,li,tf);
      k4i=dt*rhsi(t+dt,ts+k3s,ti+k3i,ta+k3a,kbar,rhoi,cpi,ls,li,ha);
      k4a=dt*rhsa(t+dt,ti+k3i,ta+k3a,rhoa,cpa,la,ha);
      ts=ts+1.0/6.0*k1s+1.0d0/3.0d0*k2s+1.0/3.0*k3s+1.0/6.0*k4s;
      ti=ti+1.0/6.0*k1i+1.0d0/3.0d0*k2i+1.0/3.0*k3i+1.0/6.0*k4i;
      ta=ta+1.0/6.0*k1a+1.0d0/3.0d0*k2a+1.0/3.0*k3a+1.0/6.0*k4a;
      t=t+dt;
      fprintf(fid,'%10.6f %10.6f %10.6f %10.6f %10.6f\r\n',t,tf,ts,ti,ta);
      outputs(j,:)=[t,tf,ts,ti,ta];
end
figure(1),clf(1)
t=outputs(:,1);
plot(t,outputs(:,2:end))