I'm trying to solve this system of ode with runge kutta 4 but the solution doesn't converge.

clearvars;close all;clc;

%parameters
q=400;
h=4;
s=2.8*3;
l=0.1;
k=1.5;
f=0.001;
c=920*850*s*l*f;
Gc=h*s;
Gd=k*s/l;
Gr=20;
Ta=20;
Ti=15;

%Define function handles
f=@(t,T) ...
         [(Gc*(Ta-T(1))+Gr*(Ta-T(1))+q*s)/c... 
         ;(Gd*(T(1)-T(2))+Gc*(Ti-T(2)))/c];
                     

%initial conditions
t(1)=1;
T(:,1)=[20,20];

%step size
h=0.001;
tfinal=100;
N=ceil(tfinal/h);
%Update loop
for i=1:N
    t(i+1)=t(i)+h;
    k1=f(t(i)    ,T(:,i));
    k2=f(t(i)+h/2,T(:,i)+h/2*k1);
    k3=f(t(i)+h/2,T(:,i)+h/2*k2);
    k4=f(t(i)+h  ,T(:,i)+h  *k3);
    T(:,i+1)=T(:,i)+h/6*(k1+2*k2+2*k3+k4);
end 
%plot the solution

figure(1); clf(1)
plot (t,T(1,:))
hold on 
plot(t,T(2,:))
xlabel('Time')
ylabel('T1 and T2')
legend('T1','T2')
grid on 
set(gcf,'color','w')

Hello,

I'm trying to solve this system of ordinary differential equations numerically by using runge kutta 4 , but the solution won't converge unless i reduce the value of parameter c by multiplying it with the factor f, even when try i reduce the time step the problem is still there, i need to solve it parameter f set to 1 , i've been struggling with this problem for day and i really appreciate if any one can help me.

[Torsten]

Did you try further reducing the stepsize h ( h=1e-7, e.g) ?

[malek chek]

Yes,i tried with h=1e-4 and it still doesn't want to converge. In the same way reducing h any further is becoming not practical.

[Torsten]

If you need f=1e-3 for convergence with h=1e-3, you'll need h=1e-6 for f=1.

I suggest you either solve your system analytically or use an implicit solver.

[malek chek]

I used ode23s and Finally, I'm getting Good results,Thank You for your help.