MATLAB: Solving a system of 6 differential equations to find a dynamic solution

differential equationsode45

I have a system of 6 different differential equations and each of these depend on each other. I have a snippet of my code as follows. I am not a pro at matlab. So my coefficients change at every time interval of 1 hour and I need to evaluate the solution at every change. However I cant seem to solve it. My output (temperature) keeps degrading incrementally.which is completely wrong. Is there any other simpler method I can follow or something wrong with my code or method? I seem to have exhausted all kinds of methods(tried manually doing RK-4 method, ode45,explicit). Any help is greatly appreciated!
My code
     z1=ah(m,1);
      z2=bh(m,1);
      z3=ch(m,1);
      z4=dh(m,1);
      z5=eh(m,1);
      z6=fh(m,1);
      z7=gh(m,1);
      z8=gne(m,1);
      z9= hh(m,1);
      z10=vh(m,1);
      z11=xh(m,1);
      z12=kh(m,1);
      z13=lh(m,1);
      z14=ph(m,1);
      z15=kill(m,1);
      z16=nec(m,1);
      z17=oh(m,1);
      z18=yh(m,1);
      z19=c22(m,1);
      MATBIT=[z1 z2 0 0 0 0;
               z4 z5 z6 0 0 0 ;
               0 z8 z9 z10 0 z11 ;
               0 0 z12 z13 z14 0 ;
               0 0 z15 z16 0 0;
               0 0 z17 0 0 z18];
      CONBIT=[z3;
               z7;
                0;
                0;
                0;
                z19];
    EQSOL1=linsolve(MATBIT,CONBIT);
    %newstate v=MATBIT*v(v=u-ue as ref to origin)
    [v,d]=eig(MATBIT);%m is value of the loop when change happens
    e1=exp(m*d(1,1));
    e2=exp(m*d(2,2));
    e3=exp(m*d(3,3));
    e4=exp(m*d(4,4));
    e5=exp(m*d(5,5));
    e6=exp(m*d(6,6));
    expmat=[e1;
            e2;
            e3;
            e4;
            e5;
            e6];
    ////////////////////////////////////////////////////////////////////
    %general solution% vt=c1*v(:,1)*expmat(1,1)+c2*v(:,2)*expmat(1,2)+c3*v(:,3)*expmat(1,3)+c4*v(:,4)*expmat(1,4)+c5*v(:,5)*expmat(1,5)+c6*v(:,6)*expmat(1,6);
    % ut=c1*v(:,1)*expmat(1,1)+c2*v(:,2)*expmat(1,2)+c3*v(:,3)*expmat(1,3)+c4*v(:,4)*expmat(1,4)+c5*v(:,5)*expmat(1,5)+c6*v(:,6)*expmat(1,6)+EQSOL1;
    % Initial value
    % um=c1*v(:,1)+c2*v(:,2)+c3*v(:,3)+c4*v(:,4)+c5*v(:,5)+c6*v(:,6)+EQSOL1;
    %/////////////////////////////////////////////////////////////////////////////
    %initial conditions
    clear um
    um=[Tinitial(m,1);
    avg1temp(1,1);
     avg1temp(1,1);
     avg1temp(1,1);
      avg1temp(1,1);
      Tin(m,1)];
    B= real(minus(um,EQSOL1));
    C=real( linsolve(v,B));
    [t,u]=ode23(@heatsystem,[0,1],um);%simulate heat through model
    tt=linspace(m-1,m,2);
    u1=real(((C(1,1)*v(1,1)*e1*tt)+(C(2,1)*v(1,2)*e2*tt)+(C(3,1)*v(1,3)*e3*tt)+(C(4,1)*v(1,4)*e4*tt)+(C(5,1)*v(1,5)*e5*tt)+(C(6,1)*v(1,6)*e6*tt)+EQSOL1(1,1)));
    u2=real(C(1,1)*v(2,1)*e1*tt+C(2,1)*v(2,2)*e2*tt+C(3,1)*v(2,3)*e3*tt+C(4,1)*v(2,4)*e4*tt+C(5,1)*v(2,5)*e5*tt+C(6,1)*v(2,6)*e6*tt+EQSOL1(2,1));
    u3=real(C(1,1)*v(3,1)*e1*tt+C(2,1)*v(3,2)*e2*tt+C(3,1)*v(3,3)*e3*tt+C(4,1)*v(3,4)*e4*tt+C(5,1)*v(3,5)*e5*tt+C(6,1)*v(3,6)*e6*tt+EQSOL1(3,1));
    u4=real(C(1,1)*v(4,1)*e1*tt+C(2,1)*v(4,2)*e2*tt+C(3,1)*v(4,3)*e3*tt+C(4,1)*v(4,4)*e4*tt+C(5,1)*v(4,5)*e5*tt+C(6,1)*v(4,6)*e6*tt+EQSOL1(4,1));
    u5=real(C(1,1)*v(5,1)*e1*tt+C(2,1)*v(5,2)*e2*tt+C(3,1)*v(5,3)*e3*tt+C(4,1)*v(5,4)*e4*tt+C(5,1)*v(5,5)*e5*tt+C(6,1)*v(5,6)*e6*tt+EQSOL1(5,1));
    u6=real(C(1,1)*v(6,1)*e1*tt+C(2,1)*v(6,2)*e2*tt+C(3,1)*v(6,3)*e3*tt+C(4,1)*v(6,4)*e4*tt+C(5,1)*v(6,5)*e5*tt+C(6,1)*v(6,6)*e6*tt+EQSOL1(6,1));
    %///////////////////////////////////////////////////////////
    %Calling each time when my coefficients values change to evaluate new solution
    function yp=heatsystem(t,y)
    global z1 z2 z3 z4 z5 z6 z7 z8 z9 z10 z11 z12 z13 z14 z15 z16 z17 z18 z19 
    yt1=z1*y(1)+z2*y(2)+z3;
    yt2=z4*y(1)+z5*y(2)+z6*y(3)+z7;
    yt3=z8*y(2)+z9*y(3)+z10*y(4)+z11*y(6);
    yt4=z12*y(3)+z13*y(4)+z14*y(5);
    yt5=z15*y(3)+z16*y(4);
    yt6=z17*y(3)+z18*y(6)+z19;
    yp=[yt1,yt2,yt3,yt4,yt5,yt6]';
    end

Best Answer

The bunch of globals called "z1, z2, ..." is really ugly. I would not dare to debug the formula. See http://www.mathworks.com/matlabcentral/answers/1971 for a reliable method to avoid globals, which are evil in general.
If you have calculated the trajectories with ode45, ode23 and a manually code RK4, and you get similar results, you can be sure, that the results are exactly, what you have coded. If you expect something else, either the code is wrong, or your expectations. You can imagine that we cannot debug e.g. the construction of MATBIT. Codes like:
      z4=dh(m,1);
      z5=eh(m,1);
      z6=fh(m,1);
      z7=gh(m,1);
or
u2=real(C(1,1)*v(2,1)*e1*tt+C(2,1)*v(2,2)*e2*tt+C(3,1)*v(2,3)*e3*tt+C(4,1)*v(2,4)*e4*tt+C(5,1)*v(2,5)*e5*tt+C(6,1)*v(2,6)*e6*tt+EQSOL1(2,1));
prevent understanding and impede debugging. So I have to say: Sorry, I do not think that anybody can find the problem. I would never trust the results of this code, even if they satisfy your expectations.
Related Question