Hi, please am trying to solve some couples of PDEs and ADEs using method of lines, but am having difficulties in converting the six matrices to 1D, PLEASE I NEED SOME , the script is attached.

% File: pde_ 1_main.m % The following equations model a plug flow reactor (PFR) % Material balance % cr(i) = ((-v*cr(i) + r1*rho)/epsilon) (1) % t z % cp(j) = ((-v*cp(j) + r2*rho)/epsilon) (2) % t z % theta_r(k) = r1- r - r3 (3) % t % theta_p(l) = r - r2 (4) % t % theta_C(m) = r3 (5) % t % theta_s(n) = r2 - r1 (6) % t % r1 =-k1*(cr(i)*(theta_s(n)) - (theta_r (k)/K1) % r =K*(theta_r(k)) % r2 =k2*(theta_p(l)) - (cp(j))*((theta_s(n))/K2) % r3 =k3*(theta_C(m))

% The variables and parameters for this model are (in cgs units) % % Reactant concentration cr(i) (eq. (1)) % Product concentration cp(j) (eq. (2)) % Surface coverage of r theta_r(k) (eq. (3)) % Surface coverage of p theta_p(l) (eq. (4)) % Coke concentration theta_C(m) (eq. (5)) % Active sites concentration theta_s(n) (eq. (6)) % Rate of adsorption reaction r1 % Rate of desorption reaction r2 % Rate of surface reaction r % Rate of deactivation reaction r3 % Rate constant for adsorption reaction k1 % Equilibrium constant for adsorption reaction K1 % Rate constant for desorption reaction k2 % Equilibrium constant for desorption reaction K2 % Rate constant for surface reaction K % Rate constant for deactivation reaction K3 % Time t % Axial position z % Entering concentration cr0 1.5 % Reactor length zl 100 % flow rate v 1 % Catalyst bulk density rho 1.0 % Void fraction epsilon 0.5

% Clear previous files clear all clc % % Global area global nz dz t cp0 theta_r0 ... z zl rho epsilon theta_p0 ... cr0 cre theta_s0 theta_C0 k1 ... v K1 K k2 K2 ... k3 cr cp nout ... theta_r theta_p theta_s theta_C y ... tf tout y0 ndzs ... ncall

% Model parameters cr0=1.5; cre=1.5; zl=100.0; v=1.0; epsilon=0.5; k1=16; K1=9; k2=30; K2=7; K=0.3; k3=2.2; cp0=0.0; theta_r0=0.0; theta_p0=0.0; theta_C0=0.0; theta_s0=2.0; rho=1.0; % Grid in axial direction nz=20; dz=zl/nz; for i=1:nz z(i)=i*dz; end

% Independent variable for ODE integration tf=200.0; tout=[0.0:50.0:tf]'; nout=5; ncall=0; % % Initial condition for i=1:nz for j=1:nz; for k=1:nz; for l=1:nz; for m=1:nz; for n=1:nz; cr(i,j,k,l,m,n)=cre; cp(i,j,k,l,m,n)=cp0; theta_r(i,j,k,l,m,n)=theta_r0; theta_p(i,j,k,l,m,n)=theta_p0; theta_C(i,j,k,l,m,n)=theta_C0; theta_s(i,j,k,l,m,n)=theta_s0; y0((i-1)*nz+(j-1)*nz*nz+(k-1)*nz*nz*nz+(l-1)*nz*nz*nz*nz+(m-1)*nz*nz*nz*nz*nz+n=cr(i,j,k,l,m,n); end end end end end end

% ODE integration reltol=1.0e-04; abstol=1.0e-04; options=odeset('RelTol',reltol,'AbsTol',abstol); mf=2; if (mf==1)[t,y]=ode15s(@pde_1,tout,y0,options);end if (mf==2)[t,y]=ode15s(@pde_2,tout,y0,options);end % % 1D to 6D matrices for it=1:nout for i=1:nz for j=1:nz; for k=1:nz; for l=1:nz; for m=1:nz; for n=1:nz; cr(i,j,k,l,m,n)=cre; cp(i,j,k,l,m,n)=cp0; theta_r(i,j,k,l,m,n)=theta_r0; theta_p(i,j,k,l,m,n)=theta_p0; theta_C(i,j,k,l,m,n)=theta_C0; theta_s(i,j,k,l,m,n)=theta_s0; cr(it,i,j,k,l,m,n)=y(it,(i-1)*nz+(j-1)*nz*nz+(k-1)*nz*nz*nz+(l-1)*nz*nz*nz*nz+(m-1)*nz*nz*nz*nz*nz+n;

end
end
end
end
end
end 
end

% % Display a heading and centerline output fprintf('\n nz = %6 d \n',nz); for it=1:nout fprintf('\n t = %4 .1f \n',t(it)); for i=1:nz fprintf('z=%5 .1f cr,cp,theta_r,theta_p,theta_C,theta_s(z,t)=%8 .5f \n',... z(i),cr,cp,theta_r,theta_p,theta_C,theta_s(it,i,1)); end end fprintf('\n ncall = %5 d \n',ncall); % % Parametric plots % % Axial profiles figure; subplot(4,1,1) plot(z,cr,cp,theta_r,theta_p,theta_C,theta_s(2,:,1));axis([0 zl 1.5 1.5]); title('cr,cp,theta_r,theta_p,theta_C,theta_s (z,t=50)'); xlabel('z'); ylabel('cr,cp,theta_r,theta_p,theta_C,theta_s (z,t=50)')

subplot(4,1,2) plot(z,cr,cp,theta_r,theta_p,theta_C,theta_s(3,:,1)); axis([0 zl 1.5 1.5]); title('cr,cp,theta_r,theta_p,theta_C,theta_s (z,t=100)'); xlabel('z'); ylabel('cr,cp,theta_r,theta_p,theta_C,theta_s (z,t=100)')

subplot(4,1,3) plot(z,cr,cp,theta_r,theta_p,theta_C,theta_s(4,:,1)); axis([0 zl 1.5 1.5]); title('cr,cp,theta_r,theta_p,theta_C,theta_s (z,t=150)'); xlabel('z'); ylabel('cr,cp,theta_r,theta_p,theta_C,theta_s (z,t=150)')

subplot(4,1,4) plot(z,cr,cp,theta_r,theta_p,theta_C,theta_s(5,:,1)); axis([0 zl 1.5 1.5]); title('cr,cp,theta_r,theta_p,theta_C,theta_s (z,t=200)'); xlabel('z'); ylabel('cr,cp,theta_r,theta_p,theta_C,theta_s (z,t=200)')