# Coupled first order differential equations

Rafique Mir
How i can solve a system of 6 first order differential equations by using numerical techniques like Euler method, RK-4th order method , ODE -45 etc.How i can solve a system of 6 first order differential equations by using numerical techniques like Euler method, RK-4th order method , ODE -45 etc.

## Answers and Replies

mathmike
can you give an example

Rafique Mir
Coupled ODEs

Dear here are the system of six differential equtions i want to solve them by Euler method , Predictor corrector method and for RK-4th order method. I have made a program in Matlab plz check wether it is right or not.
function [GG ] = myfun_rafiq(t,x)

% ----- State Variables Selection
% x1 = nw, x2 = Nw, x3 = nc, x4 = Nc, x5 = np, x6 = Np

% ----- Parameters

sigmaEjQj = 793.7;
EiQi =700;
EpQp = 13.7;
EcQc = 80;
EfQf = 0;
Vc = 9.08e6;
Vw = 1.37e7;
Vp = 1.37e6;
lk = 0.0;
Kc = 40.0;
Kp = 6.90;
f_t = 1;
sigma = 13.4e-24;
fn = 1;
fs = 0.5;
phi0 = 9.2e13;
phiE = 0.0026*phi0;
lambda = 7.4612e-005; % PER SECOND

N0 = 6.023e23;
A = 56;
S = 1.01e8;

% ---- parameters + definitio of C(t)

a = 50*3600; % in seconds
deltaCDeltaT = 20e-12;
% deltaCDeltaT=deltaCDeltaT/(3600^2);
C0 = 2.4e-13;
Cs = 25e-6;
b = 400*3600; % in seconds
t0 =50*3600; % in seconds %REPLCED 500 with 50
% t = 900*3600;
if t < a
C_t = 0;
elseif (t > a & t < b)
C_t = deltaCDeltaT*(t-t0);
else
C_t = Cs;
end
% t
% C_t

% ---- Parameters + definition of g(t)
%
tin = 500*3600;
tmax = 600*3600;
w0 = 18.3e6;
w2 = 0.3*w0;
alpha = 0.004;
% g_t = 1;

if t <= tin
g_t = 1;
elseif (t > tin & t <= tmax)
g_t = 1 - alpha*(tin - t);
else %REPLACED 'b' with tmax
g_t = w2/w0;
end
% if(t<tin)
% g_t = 1;
% elseif((t<=tmax))
% g_t = 1 - alpha*(tin - t);
% % g_t=0.00001;
% else %REPLACED 'b' with tmax
% g_t = w2/w0;
% % g_t=0.00001;
% end

% % pause

Sw = (C_t* S*N0*fn*fs)/(Vw*A);

dxdt = zeros(6,1);

dxdt(1) = sigma*phiE*x(2) - ( (sigmaEjQj*g_t/Vw) + lk*g_t/Vw + lambda)*x(1) + (Kp*g_t/Vw)*x(5) + (Kc*g_t/Vw)*x(3);

dxdt(2) = -( (sigmaEjQj*g_t/Vw) + lk*g_t/Vw + sigma*phiE)*x(2) + (Kp*g_t/Vw)*x(6) + (Kc*g_t/Vw)*x(4) + Sw;

dxdt(3) = sigma*phi0*x(4) + (EcQc*g_t/Vc)*x(1) - ( Kc*g_t/Vc + lambda)*x(3);

dxdt(4) = (EcQc*g_t/Vc)*x(2) - ( Kc*g_t/Vc + sigma*phi0)*x(4);

dxdt(5) = (EpQp*g_t/Vp)*x(1) - ( Kp*g_t/Vp + lambda)*x(5);

dxdt(6) = (EpQp*g_t/Vp)*x(2) - ( Kp*g_t/Vp)*x(6);

GG=[dxdt(1); dxdt(2); dxdt(3); dxdt(4); dxdt(5); dxdt(6)];

% ---- PIEAS

% Clearing work space

clc;
clear all;
% close all;
tic
% Define time for simulation
%x0 = [0.1, 0.2, 0.3, 0.4, 0.1, 0.1]*10^-7;
% x0=[0 0.5 0 0.5 0 0.5];
% T = 1:500*3600; % 10 seconds
% T=T*3600;

% Define initional conditions

%R=myfun_rahat(T,x0);

% Run simulation
% options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);

[tt,xx] = ode45(@myfun_rafiq,[0 1000*3600],[0 0.5 0 0.5 0 50]);

toc

% Plot Results
grid on;
figure(1),plot( tt, xx(:,1),'r.-',tt,xx(:,3),'m:',tt,xx(:,5),'b.-');
xlabel('t(sec)');
title('Specific Activity');
legend('nw','nc','np');
figure(2)
grid on;
plot(tt, xx(:,1),'r-');
%axis([0 1.81e5 0 3e-4]);

Homework Helper
Rafique Mir said:
Dear here are the system of six differential equtions i want to solve them by Euler method , Predictor corrector method and for RK-4th order method. I have made a program in Matlab plz check wether it is right or not.

dxdt(1) = sigma*phiE*x(2) - ( (sigmaEjQj*g_t/Vw) + lk*g_t/Vw + lambda)*x(1) + (Kp*g_t/Vw)*x(5) + (Kc*g_t/Vw)*x(3);

Oh my God dude! Now, I really admire you for taking the time to input all that code but it's just too hard to read. I picked out the first one above and I can't understand it. I realize you're probably not familiar with the math-formatting code we use in here called "LaTex" but that would help with understanding your code. If you want to try using it, you can go to the General Physics forum and read the "Introducing LaTex" thread at the top. Also, I use Mathematica so would not be able to help you with Matlab.