Problem in Iteration or model

Hi,

I have written a program to solve bunch of ODEs and during dynamic run (space and time), the output values do not change with time. Output only changes with space. Any help is appreciable. Thanks in advance.

My program is as follows:

%                                   In The Name Of ALLAH swt
% ---------------------------------------------------------------------------------------------
% This m-file is the main one for solving the set of PDEs of the model.
% For solving the set of PDEs, the equations have been discretized in z and sloved dynamically
% in each point along the reactor.

clear all close all SteadySMSR clc global Pt0 u delz deltaP u01 u02 u03 u04 u05 u06 u013 u015 u016 u0 counter s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 kg hf r1 r2 r3

u0=[2.04e-1;4.513e-2;1.33e-2;6.682e-1;1.312e-1;9.75e-3;2.06e-1;4.58e-2;1.36e-2;6.72e-1;1.33e-1;9.75e-3;930;935;4.082;1];

u01=2.50e-1; u02=7.47e-3; u03=2.22e-3; u04=7.22e-1; u05=8.34e-3; u06=9.75e-3; u013=893; u015=4.083; %***Inlet condition*****

Pt0=24; L=14.35; % m no_incrmnt=input('Please specify the no. of incriments:\n'); delz=L/no_incrmnt; % m

u=u0; tfinal=1200; tinc=10; tmesh=1:tinc:tfinal;

s1(1,1:length(tmesh))=2.50e-1; s2(1,1:length(tmesh))=7.47e-3; s3(1,1:length(tmesh))=2.22e-3; s4(1,1:length(tmesh))=7.22e-1; s5(1,1:length(tmesh))=8.34e-3; s6(1,1:length(tmesh))=9.75e-3; s7(1,1:length(tmesh))=2.50e-1; s8(1,1:length(tmesh))=7.47e-3; s9(1,1:length(tmesh))=2.22e-3; s10(1,1:length(tmesh))=7.22e-1; s11(1,1:length(tmesh))=8.34e-3; s12(1,1:length(tmesh))=9.75e-3; s13(1,1:length(tmesh))=893; s14(1,1:length(tmesh))=893; s15(1,1:length(tmesh))=4.083; s16(1,1:length(tmesh))=1;

f1='Report1'; f2='Report2'; f3='Report3'; f4='Report4'; f5='Report5'; f6='Report6'; f7='Report7'; f8='Report8'; f9='Report9'; f10='Report10'; f11='Report11'; f12='Report12'; f13='Report13'; f14='Report14'; f15='Report15'; f16='Report16';

options=optimset('LargeScale','off','TolX',1e-10,'MaxIter',1e9,'MaxFunEvals',2e20);

counter=2; waithandle=waitbar(0,' Step2 : ------ Dynamic Simulation ------ '); % Initialize waitbar

% ----- Solving the system of PDEs -----

while counter<(no_incrmnt+2) if (counter+0.01<=(no_incrmnt+2)), waitbar(counter/(no_incrmnt+1),waithandle); end

    [t1 u1]=ode23s(@simulationMSR,1:tinc:tfinal-1,u0,options);

    s1(counter,:)=u1(:,1)';
    s2(counter,:)=u1(:,2)';
    s3(counter,:)=u1(:,3)';
    s4(counter,:)=u1(:,4)';
    s5(counter,:)=u1(:,5)';
    s6(counter,:)=u1(:,6)';
    s7(counter,:)=u1(:,7)';
    s8(counter,:)=u1(:,8)';
    s9(counter,:)=u1(:,9)';
    s10(counter,:)=u1(:,10)';
    s11(counter,:)=u1(:,11)';
    s12(counter,:)=u1(:,12)';
    s13(counter,:)=u1(:,13)';
    s14(counter,:)=u1(:,14)';
    s15(counter,:)=u1(:,15)';
    s16(counter,:)=u1(:,16)';

    load('ReportS','w');

    a1=w(counter,1);
    a2=w(counter,2);
    a3=w(counter,3);
    a4=w(counter,4);
    a5=w(counter,5);
    a6=w(counter,6);
    a7=w(counter,7);
    a8=w(counter,8);
    a9=w(counter,9);
    a10=w(counter,10);
    a11=w(counter,11);
    a12=w(counter,12);
    a13=w(counter,13);
    a14=w(counter,14);
    a15=w(counter,15);
    a16=1;

    u0=[a1;a2;a3;a4;a5;a6;a7;a8;a9;a10;a11;a12;a13;a14;a15;a16];
    R1(counter-1)=r1; R2(counter-1)=r2; R3(counter-1)=r3;
    Pt0=Pt0-(deltaP/L)*delz;
    counter=counter+1;

end close(waithandle);

% ----- Saves the results -----

save(f1,'s1'); save(f2,'s2'); save(f3,'s3'); save(f4,'s4'); save(f5,'s5'); save(f6,'s6'); save(f7,'s7'); save(f8,'s8'); save(f9,'s9'); save(f10,'s10'); save(f11,'s11'); save(f12,'s12'); save(f13,'s13'); save(f14,'s14'); save(f15,'s15'); save(f16,'s16');

bbb=clock;

% ----- 3D plotts of the model -----

figure(1); x=1:1:(no_incrmnt+1); time=1:tinc:tfinal; z=0:delz:L;

figure(1);mesh(time,z,s1);xlabel('Time (days)'); figure(2);mesh(time,z,s2);xlabel('Time (days)'); figure(3);mesh(time,z,s3);xlabel('Time (days)'); figure(4);mesh(time,z,s4);xlabel('Time (days)'); figure(5);mesh(time,z,s5);xlabel('Time (days)'); % figure(6);mesh(time,x,s6); figure(7);mesh(time,z,s13); figure(8);mesh(time,z,s15); s16(1,:)=s16(2,:) figure(9);mesh(time,z,s16);