# MATLAB: Runge kutta 4th order

ode questionrunge kuttarunge kutta 4th order

function Untitled% Runge kutta 4th order method% dCa/dt = F(Ca,in - Ca)/V -2k(T)Ca^2-----eq.1% dT/dt = F(Tin - T)/V +(-deltaH/rhoCp)*2*k(T)Ca^2-(T-Tf)UA/VrhoCp----eq.2clc; clear;%constantsFlowrate = 20;V = 100;U_A = 20000;Cp = 4.2;rho = 1000;delta_H = 596619;Ar = 6.85E+11;E = 76534.704;R = 8.314;T_j = 250;%initial conditions t_initial(1) = 0;T_initial(1) = 275;Ca_initial(1) = 1;%Arrhenius equationn = 10 ;   for j = 1:n   k(j) = Ar*exp(-E/(R*T_initial(j)));   if j<n   T_initial(j+1) = T_initial(j) +5;   endend%Define Function handlefCa = @(t,Ca,T) Flowrate*(((Ca_initial - Ca)/(V)) -(2*k(j)*Ca*Ca));fT  = @(t,Ca,T) Flowrate*(T_initial- T)/V +(-delta_H/(rho*Cp))*2*k(j)*Ca*Ca-(T-T_j)*(U_A/(V*rho*Cp));%step sizeh = 10;t_f = 300; % final timeN = (t_f/h);% update loopfor i = 1:N    t_initial(i+1) = t_initial(i) + h;    K1Ca_initial = fCa (t_initial(i)      ,Ca_initial(i)                   ,T_initial(i)                 );    K1T_initial  = fT  (t_initial(i)      ,Ca_initial(i)                   ,T_initial(i)                 );    K2Ca_initial = fCa (t_initial(i)+ h/2 ,Ca_initial(i)+ h/2*K1Ca_initial ,T_initial(i)+ h/2*K1T_initial);    K2T_initial  = fT  (t_initial(i)+ h/2 ,Ca_initial(i)+ h/2*K1Ca_initial ,T_initial(i)+ h/2*K1T_initial);    K3Ca_initial = fCa (t_initial(i)+ h/2 ,Ca_initial(i)+ h/2*K2Ca_initial ,T_initial(i)+ h/2*K2T_initial);    K3T_initial  = fT  (t_initial(i)+ h/2 ,Ca_initial(i)+ h/2*K2Ca_initial ,T_initial(i)+ h/2*K2T_initial);    K4Ca_initial = fCa (t_initial(i)+ h   ,Ca_initial(i)+ h  *K3Ca_initial ,T_initial(i)+ h  *K3T_initial);    K4T_initial  = fT  (t_initial(i)+ h   ,Ca_initial(i)+ h  *K3Ca_initial ,T_initial(i)+ h  *K3T_initial);        Ca_initial(i+1) = Ca_initial(i) +h/6*(K1Ca_initial + 2*K2Ca_initial +  2*K3Ca_initial + K4Ca_initial);    T_initial (i+1) = T_initial(i)  +h/6*(K1T_initial  + 2*K2T_initial  +  2*K3T_initial  + K4T_initial );endplot(t_initial,Ca_initial)hold onplot(t_initial,T_initial)xlabel('time')ylabel('Concentration')legend ('Conc.' , 'Temp.')set(gca,'fontsize',16)
———————————————————————
I'm getting error 'T_initial (i+1) = T_initial(i) +h/6*(K1T_initial + 2*K2T_initial + 2*K3T_initial + K4T_initial );' here. It's saying 'nable to perform assignment because the left and right sides have a different number of elements.' where am i going wrong ? Please help
Here k(T) is k is function of T i.e temperature and k is of the arrhenius equation.

You are using k(j) in your derivative function handles. This is a mistake since it causes the function handles to use the current value of the k(j) element at the time of function handle creation as a constant in the function handle. I.e., whatever the last value of k(j) you calculated in the loop filling up the k vector will be used for every derivative call.
Another mistake is making T_initial a vector and using it in a function handle, creating a vector result. This will cause size mismatches in assignments, etc.
You have two variables, Ca and T, so you need two scalar derivative functions for these. And since k is a function of T, you need to calculate k on the fly using the current value of T ... you don't calculate this ahead of time since you don't know T ahead of time.
It is not entirely clear from what you have written, but I would say that these equations
% dCa/dt = F(Ca,in - Ca)/V -2k(T)Ca^2-----eq.1
% dT/dt = F(Tin - T)/V +(-deltaH/rhoCp)*2*k(T)Ca^2-(T-Tf)UA/VrhoCp----eq.2
would translate into code like this
k = @(T)Ar*exp(-E/(R*T)); % not sure about thisfCa = @(t,Ca,T) Flowrate*(((Ca_initial - Ca)/(V)) -(2*k(T)*Ca*Ca));fT  = @(t,Ca,T) Flowrate*(T_initial- T)/V +(-delta_H/(rho*Cp))*2*k(T)*Ca*Ca-(T-T_j)*(U_A/(V*rho*Cp));
and you would get rid of this code entirely:
n = 10 ;   % delete thisfor j = 1:n   % delete this   k(j) = Ar*exp(-E/(R*T_initial(j)));   % delete this   if j<n   % delete this   T_initial(j+1) = T_initial(j) +5;   % delete this   end   % delete thisend   % delete this
But some of this is guesswork on my part, particularly the k(T) stuff, because you haven't posted enough information. Can you post an image of the actual differential equations you are trying to solve so we can verify this?
I would also strongly advise that you put comments on every line with a constant and document the units of that constant. This makes it much easier to debug and spot unit errors in the code, and much easier for reviewers to know what is going on in your code.