System of mixed PDE with 2 variables using FDM and RK4

I want solve the attached mixed pdes in 2 variables with FDM and R4K. See my trial below: I received and error that 'N' is uncrecognised.

% parameters.

k1 = 6*10^(-12);

eta1 = 0.0240;

apha3 = -0.001104;

gama1 = 0.1093;

d = 0.0002;

%thetab = 0.0001;

% activity parameter

xi = 0.1;

%%%% Simplify parameters

A = (apha3*k1)/gama1;

B = A/k1;

%%%% Discretize xspace

z_vec = linspace(0, d, N);

dz = z_vec(2) - z_vec(1);

%%%% Discretize t

dt = (dz^2)/(2);

%%%% Discretize t

dtrk4 = (dz^2)/(2);

trk4_vec = 0: dtrk4:10;

% Allocate memory for u

thetark4_mat = zeros(length(z_vec), length(trk4_vec));

vrk4_mat = zeros(length(z_vec), length(trk4_vec));

%%%% IC

thetark4_mat(1,:) = 0; % Left end of pipe

thetark4_mat(end,:) = d; % right end of the pipe

vrk4_mat(1,:) = 0; % Left end of pipe

vrk4_mat(end,:) = d; % right end of the pipe

for tdz = 1:length(trk4_vec)-1

k11 = Derivative(thetark4_mat(:, tdz));

k12 = Derivative(vrk4_mat(:, tdz));

k21 = Derivative(thetark4_mat(:, tdz) + k11*dtrk4/2);

k22 = Derivative(vrk4_mat(:, tdz) + k12*dtrk4/2);

k31 = Derivative(thetark4_mat(:, tdz) + k21*dtrk4/2);

k32 = Derivative(vrk4_mat(:, tdz) + k22*dtrk4/2);

k41 = Derivative(thetark4_mat(:, tdz) + k31*dtrk4);

k42 = Derivative(vrk4_mat(:, tdz) + k32*dtrk4);

phi = (1/6)*(k11+ 2*k21 + 2*k31 +k41);

phi1 = (1/6)*(k12+ 2*k22 + 2*k32 +k42);

thetark4_mat(:, tdz+1) = thetark4_mat(:, tdz) + phi*dtrk4;

vrk4_mat(:, tdz+1) = vrk4_mat(:, tdz) + phi*dtrk4;

end

%%%% Plot this

[tt, zz] = meshgrid(trk4_vec, z_vec);

mesh(zz, tt, thetark4_mat)

xlabel('z')

ylabel('time')

zlabel('theta')

figure()

mesh(zz, tt, vrk4_mat)

xlabel('z')

ylabel('time')

zlabel('v')

function dthetadt = Derivative(theta, v)

%global x_vec dx k

dthetadt = 0*theta;

dvdt = 0*theta;

z_vec = linspace(0,10,10);

% parameters.

k1 = 6*10^(-12);

eta1 = 0.0240;

apha3 = -0.001104;

gama1 = 0.1093;

d = 0.0002;

%thetab = 0.0001;

% activity parameter

xi = 0.1;

%%%% Simplify parameters

A = (apha3*k1)/gama1;

B = A/k1;

dz = z_vec(2) - z_vec(1);

for j = 2:length(z_vec)-1

dvdt(j) = 0;

dthetadt(j) = A/(apha3*dz^2)*(theta(j+1) -2*theta(j) + theta(j-1)) + ...

B/(2*dz)*(v(j+1) - v(j-1));

dvdt(j) = A/(2*dz^3)*(theta(j+2) +2*theta(j+1) -2*theta(j-1)- theta(j-2))+ ...

(eta1 - B)/(dz^2)*(v(j+1) -2*v(j) + v(j-1)) - z/(2*dz)*(theta(j+1) - theta(j-1));

end

end

[Torsten]

You try to solve

dv/dt = d/dz(eta1*dv/dz + alpha3*K1/gamma1*d^2theta/dz^2 - alpha3/gamma1*dv/dz + zeta*theta)

not

0 = d/dz(eta1*dv/dz + alpha3*K1/gamma1*d^2theta/dz^2 - alpha3/gamma1*dv/dz + zeta*theta)

And why do you use your own code and not one from the MATLAB suite ?

And the error message is clear:

You didn't define N before you used it in the command

z_vec = linspace(0, d, N);

[University Glasgow]

I couldn't get the pde or ode solver that could solve my kind of system. That's why I tried to use FDM and RK4. But I will glad if you could provide another code. Thank you.

[University Glasgow]

From the second equation, the dv/dt component is zero.

[Torsten]

Yes, it should be zero, but in your code, you don't solve it as to be zero.

At the moment, you have a usual system of two discretized pdes that can easily be transfered to ODE15S.

What specialities do you plan that make it only solvable with your individual RK4 code ?

As I said: Make a literature research concerning a suitable discretization scheme - I don't know either which way to go best for these equations.

[University Glasgow]

Thank you for suggesting ODE15S.

[University Glasgow]

Hi, I have been able to use the ode15s however, I still receive some error. For a start, I used N = 5. Please, how do I write the code for the general case. In the attached file, h = d/N. I discretise by setting z_i = ih. The initial conditions are all 0s.

[Torsten]

But you already wrote the function "Derivative" for the general case in your RK4 code, didn't you ?

[University Glasgow]

It didn’t work. Please help me out

[Torsten]

I think from here on it's again your turn.

u0 = [0; 0; 0; 0; 0];
v0 = [0; 0; 0; 0; 0];
tspan = [0 0.0002];
M = [0 0 0 0 0 0 0 0 0 0; 0 1 0 0 0 0 0 0 0 0; 0 0 1 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 ; 0 0 0 0 0 0 0 0 0 0;...
    0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0];
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
[t,y] = ode15s(@lcode1,tspan,[u0;v0], options);

N = 5

Index exceeds the number of array elements. Index must not exceed 5.

Error in solution>lcode1 (line 40)
(G/(h^3))*(-u(3) + 3*u(4)-3*u(5)-u(6))+(C/(2*h^2))*(v(4) -2*v(3) - v(2)) + (xi/(2*h))*(u(4)-u(2))

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode15s (line 153)
odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

u = y(:,1:5)
v = y(:,6:10)
%%
% Plot the solution.

function rhsode = lcode1(t, y)
% parameters.
k1 = 6*10^(-12);    
eta1 = 0.0240;
apha3 = -0.001104;
gama1 = 0.1093;
d = 0.0002;
N =5
% activity parameter
xi = 0.1; 
h = d/N;   %% step size

%%%% Simplify parameters
A = k1/gama1;
B = apha3/k1;
C = eta1 - B;
G = apha3*A;
Phi = 0;
u = y(1:5);
v = y(6:10);
rhsode = [u(1) -Phi
    (A/(h^2))*(u(3)-2*u(2)-u(1))+(B/(2*h))*(v(3) -v(1))
    (A/(h^2))*(u(4)-2*u(3)-u(2))+(B/(2*h))*(v(4) -v(2))
    (A/(h^2))*(u(5)-2*u(4)-u(3))+(B/(2*h))*(v(5) -v(3))
    u(5) -Phi
    v(1)
    (G/(h^3))*(-u(2) + 3*u(3)-3*u(4)-u(5))+(C/(2*h^2))*(v(3) -2*v(2) - v(1)) + (xi/(2*h))*(u(3)-u(1))
    (G/(h^3))*(-u(3) + 3*u(4)-3*u(5)-u(6))+(C/(2*h^2))*(v(4) -2*v(3) - v(2)) + (xi/(2*h))*(u(4)-u(2))
    (G/(h^3))*(-u(4) + 3*u(5)-3*u(6)-u(7))+(C/(2*h^2))*(v(5) -2*v(4) - v(3)) + (xi/(2*h))*(u(5)-v(3))
    v(5)];
end

[University Glasgow]

Okay, let me try thank you.

[University Glasgow]

Find attached the refined, code but I still received this error: Error in lcmodel (line 6) [t,y] = ode15s(@lcode1,tspan,[u0;v0], options);

I tried to fixed the: "Index exceeds the number of array elements. Index must not exceed 5" by setting N = 4.

Changes made: All I did was to created a function matrixode that returns the matrix M and I changed N from 5 to 4.

[Torsten]

I thought the error message after my changes was clear: You address u(6) and u(7), but the vector u has only 5 elements:

(G/(h^3))*(-u(3) + 3*u(4)-3*u(5)-u(6))+(C/(2*h^2))*(v(4) -2*v(3) - v(2)) + (xi/(2*h))*(u(4)-u(2))

(G/(h^3))*(-u(4) + 3*u(5)-3*u(6)-u(7))+(C/(2*h^2))*(v(5) -2*v(4) - v(3)) + (xi/(2*h))*(u(5)-v(3))

[University Glasgow]

Ah, I see, I think is from the choice of the series of the d^3theta/dz^3 discretization. Thank you will fix it and try it.

[University Glasgow]

Despite the deleting the u(6) and u(7) I still receive this error: Error in lcmodel (line 6) t,y] = ode15s(@lcode1,tspan,[u0;v0], options);

[Torsten]

This error message means that ODE15S has problems with the solution of your algebraic equations for v.

I gave you the advice to make a literature research for the discretization because I foresaw that this error would appear using a blind-eye approach.

I can't help you from here any further before the problem of discretization is solved satisfactorily.

u0 = [0; 0; 0; 0; 0];
v0 = [0; 0; 0; 0; 0];
tspan = [0 0.0002];
M =  [0 0 0 0 0 0 0 0 0 0; 0 1 0 0 0 0 0 0 0 0; 0 0 1 0 0 0 0 0 0 0; 0 0 1 0 0 0 0 0 0 0 ; 0 0 0 0 0 0 0 0 0 0;...
    0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0];
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
[t,y] = ode15s(@lcode1,tspan,[u0;v0], options);

Error using daeic12
This DAE appears to be of index greater than 1.

Error in ode15s (line 305)
[y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(ode,odeArgs,t,ICtype,Mt,y,yp0,f0,...

u = y(:,1:5)
v = y(:,6:10)
% Plot the solution.
function rhsode = lcode1(t,y)
% parameters.
k1 = 6*10^(-12);    
eta1 = 0.0240;
apha3 = -0.001104;
gama1 = 0.1093;
d = 0.0002;
N = 5;
% activity parameter
xi = 0.1; 
h = d/N;   %% step size
%%%% Simplify parameters
A = k1/gama1;
B = apha3/k1;
C = eta1 - B;
G = apha3*A;
Phi = 0;
u = y(1:5);
v = y(6:10);
rhsode = [u(1)-Phi
         (A/(h^2))*(u(3)-2*u(2)-u(1))+(B/(2*h))*(v(3)-v(1))
         (A/(h^2))*(u(4)-2*u(3)-u(2))+(B/(2*h))*(v(4)-v(2))
         (A/(h^2))*(u(5)-2*u(4)-u(3))+(B/(2*h))*(v(5)-v(3))
         u(5)-Phi
         v(1)
         (G/(h^3))*(-u(2)+3*u(3)-3*u(4))+(C/(2*h^2))*(v(3)-2*v(2)-v(1))+(xi/(2*h))*(u(3)-u(1))
         (G/(h^3))*(-u(3)+3*u(4)-3*u(5))+(C/(2*h^2))*(v(4)-2*v(3)-v(2))+(xi/(2*h))*(u(4)-u(2))
         (G/(h^3))*(-u(4)+3*u(5))+(C/(2*h^2))*(v(5)-2*v(4)-v(3))+(xi/(2*h))*(u(5)-v(3))
         v(5)];
end

[Torsten]

Here is a way to circumvent the problem with the algebraic equations for v:

u0 = [0; 0; 0; 0; 0];
%v0 = [0; 0; 0; 0; 0];
tspan = [0 0.0002];
%M =  [0 0 0 0 0 0 0 0 0 0; 0 1 0 0 0 0 0 0 0 0; 0 0 1 0 0 0 0 0 0 0; 0 0 1 0 0 0 0 0 0 0 ; 0 0 0 0 0 0 0 0 0 0;...
    %0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0];
M = [0 0 0 0 0;0 1 0 0 0;0 0 1 0 0;0 0 0 1 0;0 0 0 0 0];
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
%[t,y] = ode15s(@lcode1,tspan,[u0;v0], options);
[t,y] = ode15s(@lcode1,tspan,u0, options);

u = y(:,1:5)

u = 11×5

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

%v = y(:,6:10)
% Plot the solution.
function rhsode = lcode1(t,y)
% parameters.
k1 = 6*10^(-12);    
eta1 = 0.0240;
apha3 = -0.001104;
gama1 = 0.1093;
d = 0.0002;
N = 5;
% activity parameter
xi = 0.1; 
h = d/N;   %% step size
%%%% Simplify parameters
A = k1/gama1;
B = apha3/k1;
C = eta1 - B;
G = apha3*A;
Phi = 0;
u = y(1:5);
%v = y(6:10);
fun = @(v)[v(1);...
           (G/(h^3))*(-u(2)+3*u(3)-3*u(4))+(C/(2*h^2))*(v(3)-2*v(2)-v(1))+(xi/(2*h))*(u(3)-u(1));...
           (G/(h^3))*(-u(3)+3*u(4)-3*u(5))+(C/(2*h^2))*(v(4)-2*v(3)-v(2))+(xi/(2*h))*(u(4)-u(2));...
           (G/(h^3))*(-u(4)+3*u(5))+(C/(2*h^2))*(v(5)-2*v(4)-v(3))+(xi/(2*h))*(u(5)-v(3));...
           v(5)];
options = optimset('Display','off');
v = fsolve(fun,zeros(5,1),options);
rhsode = [u(1)-Phi
         (A/(h^2))*(u(3)-2*u(2)-u(1))+(B/(2*h))*(v(3)-v(1))
         (A/(h^2))*(u(4)-2*u(3)-u(2))+(B/(2*h))*(v(4)-v(2))
         (A/(h^2))*(u(5)-2*u(4)-u(3))+(B/(2*h))*(v(5)-v(3))
         u(5)-Phi];
         %v(1)
         %(G/(h^3))*(-u(2)+3*u(3)-3*u(4))+(C/(2*h^2))*(v(3)-2*v(2)-v(1))+(xi/(2*h))*(u(3)-u(1))
         %(G/(h^3))*(-u(3)+3*u(4)-3*u(5))+(C/(2*h^2))*(v(4)-2*v(3)-v(2))+(xi/(2*h))*(u(4)-u(2))
         %(G/(h^3))*(-u(4)+3*u(5))+(C/(2*h^2))*(v(5)-2*v(4)-v(3))+(xi/(2*h))*(u(5)-v(3))
         %v(5)];
end

[University Glasgow]

Thank you very much.

[University Glasgow]

Please can you check for me? Error from the ode solver.

[Torsten]

theta and v both have only 5 elements.

But in these loops

for i=(N+4):(2*N)
    rhsode(i,1)= C/(2*h^3)*(v(i+1) + 2*v(i) - 2*v(i-1))+ ...
        + G/(h^3)*(theta(i+2) -2*theta(i+1) + 2*theta(i-1) - theta(i-2))...
        +xi/(2*h)*(theta(i+1) - theta(i-1)); %uses central difference for third derivative for theta
end

rhsode(2*N+1,1)=0;
for i = (5+N):(2*N + 1)
    rhsode(i,1)= C/(2*h^3)*(v(i+1) + 2*v(i) - 2*v(i-1))+ ...
        + G/(h^3)*(-theta(i-2) +3*theta(i-1) -3*theta(i) + theta(i+1))...
        +xi/(2*h)*(theta(i+1) - theta(i-1)); %involves special third derivative
end

almost every index of theta and v exceeds 5.

Use two arrays rhsode_theta and rhsode_v both arrays having 5 elements. Then indexing should be simple.

At the end, return

rhsode = [rhsode_theta;rhsode_v]

[University Glasgow]

Oh, thank you

Answer 1

[University Glasgow]

% I still don't this error despite having 5 by 5 entries in both you theta and v:

Error in counts_4_IjuptilK_080822 (line 41)

[t,y] = ode15s(@lcode1,tspan,u0, options);

% parameters.

k1 = 6*10^(-12);

eta1 = 0.0240;

apha3 = -0.001104;

gama1 = 0.1093;

d = 0.0002;

N = 4;

% Boundary Conditions

Phi = 0;

% activity parameter

xi = 0.1;

h = d/N; %% step size

%%%% Simplify parameters

A = k1/gama1;

B = apha3/k1;

C = eta1 - B;

G = apha3*A;

Theta = 0.0001;

z=linspace(0,d,N+1);

theta0 = Theta*sin(pi*z/d)

v0 = zeros(1,N+1)

M1=eye(N+1,N+1);

M1(1,1)=0;

M1(N+1,N+1)=0;

M2=zeros(N+1,N+1);

M=[M1 M2;M2 M2];

u0 = [theta0'; v0'];

tspan = [0 10];

options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);

%[t,y] = ode15s(@lcode1,tspan,[u0;v0], options);

[t,y] = ode15s(@lcode1,tspan,u0, options);

%theta = y(:,1:5);

%v = y(:,1:5);

% Plot the solution.

function rhsode = lcode1(~,y, k1, eta1, apha3, d, Phi, xi, A, B, C, G)

theta = y(1:5);

v = y(6:10);

rhsode = [theta(1)-Phi

(A/(h^2))*(theta(3)+2*theta(2)-theta(1))+(B/(2*h))*(v(3)-v(1))

(A/(h^2))*(theta(4)+2*theta(3)-theta(2))+(B/(2*h))*(v(4)-v(2))

(A/(h^2))*(theta(5)+2*theta(4)-theta(3))+(B/(2*h))*(v(5)-v(3))

theta(5)-Phi

v(N+1)

(G/(h^3))*(-theta(5) + 3*theta(2)-3*theta(3)+theta(4))+(C/(2*h^2))*(v(3) +2*v(2) - v(1)) + (xi/(2*h))*(theta(3)-theta(1))

(G/(h^3))*(-theta(2) + 3*theta(3)-3*theta(4)+theta(5))+(C/(2*h^2))*(v(4) -2*v(3) - v(2)) + (xi/(2*h))*(theta(4)-theta(2))

(G/(h^3))*(-theta(3) + 3*theta(4)-3*theta(5)+theta(2))+(C/(2*h^2))*(v(5) -2*v(4) - v(3)) + (xi/(2*h))*(theta(5)-theta(3))

v(5)];

end

[Torsten]

I don't understand why you deleted the call to "fsolve" in an earlier code.

Here is the error again I already commented on.

% parameters.
k1 = 6*10^(-12);     
eta1 = 0.0240;
apha3 = -0.001104;
gama1 = 0.1093;
d = 0.0002;
N = 4;

% Boundary Conditions
Phi = 0;
% activity parameter
xi = 0.1; 
h = d/N;   %% step size
%%%% Simplify parameters
A = k1/gama1;
B = apha3/k1;
C = eta1 - B;
G = apha3*A;
Theta = 0.0001;
z=linspace(0,d,N+1);
theta0 = Theta*sin(pi*z/d)
v0 = zeros(1,N+1)
M1=eye(N+1,N+1);
M1(1,1)=0;
M1(N+1,N+1)=0;
M2=zeros(N+1,N+1);
M=[M1 M2;M2 M2];
u0 = [theta0'; v0'];
tspan = [0 10];
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
%[t,y] = ode15s(@lcode1,tspan,[u0;v0], options);
[t,y] = ode15s(@(t,y)lcode1(t,y,k1, eta1, apha3, d, Phi, h, xi, A, B, C, G),tspan,u0, options);
%theta = y(:,1:5);
%v = y(:,1:5);
% Plot the solution.
function rhsode = lcode1(~,y, k1, eta1, apha3, d, Phi, h, xi, A, B, C, G)
theta = y(1:5);
v = y(6:10);
rhsode = [theta(1)-Phi 
    (A/(h^2))*(theta(3)+2*theta(2)-theta(1))+(B/(2*h))*(v(3)-v(1))
    (A/(h^2))*(theta(4)+2*theta(3)-theta(2))+(B/(2*h))*(v(4)-v(2))
    (A/(h^2))*(theta(5)+2*theta(4)-theta(3))+(B/(2*h))*(v(5)-v(3))
    theta(5)-Phi 
    v(1)
    (G/(h^3))*(-theta(5) + 3*theta(2)-3*theta(3)+theta(4))+(C/(2*h^2))*(v(3) +2*v(2) - v(1)) + (xi/(2*h))*(theta(3)-theta(1))
    (G/(h^3))*(-theta(2) + 3*theta(3)-3*theta(4)+theta(5))+(C/(2*h^2))*(v(4) -2*v(3) - v(2)) + (xi/(2*h))*(theta(4)-theta(2))
    (G/(h^3))*(-theta(3) + 3*theta(4)-3*theta(5)+theta(2))+(C/(2*h^2))*(v(5) -2*v(4) - v(3)) + (xi/(2*h))*(theta(5)-theta(3))         
    v(5)];
end

Answer 2

[University Glasgow]

Hi, thank you for your help. The code is working well for 4 counts.

But when I tried to generalised, I received this error: Error using daeic12 This DAE appears to be of index greater than 1.

Below is my code:

clear

clc

close all

% parameters.

global k1 eta1 alpha3 gamma1 d N Phi xi h A B C G

k1 = 6*10^(-12);

eta1 = 0.0240;

xi = -0.1; % activity parameter

alpha3 = -0.001104;

gamma1 = 0.1093;

d = 0.0002;

N = 4;

% Boundary Conditions

Phi = 0;

%%%% Simplify parameters

A = k1/gamma1;

B = alpha3/k1;

C = eta1 - B;

G = alpha3*A;

% step size

h = d/N;

% range of z

z=linspace(0,d,N+1);

% initials

Theta = 0.0001;

theta0 = Theta*sin(pi*z/d);

v0 = zeros(1,N+1);

theta0_int=theta0(1:N);

v0_int=v0(1:N);

u0 = [theta0_int'; v0_int'];

% Matrix M

M1=eye(N,N);

M2=zeros(N,N);

M=[M1 M2;M2 M2];

% t span

tspan = [0 4];

% ode solver

options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);

[t,y] = ode15s(@lcode1,tspan,u0, options);

% Extract the solution for theta and v

%theta = [Phi*ones(length(t), 1) y(:,1:N-1) Phi*ones(length(t), 1)];

%v = [zeros(length(t), 1) y(:,N:2*(N-1)) zeros(length(t), 1)];

% Right hand side of the ODEs: F(U)

function rhsode = lcode1(t,y)

global Phi xi h A B C G N

% initialize theta and v

theta = y(1:N);

v = y(N+1:2*N);

% theta equation

% discretise equation:(A/(h^2))*(theta(i+1)+2*theta(i)-theta(i-1))+(B/(2*h))*(v(i+1)-v(i-1))

for i=2:N

theta(N+1) = 0;

v(N+1) = 0;

rhsode(i,1)=(A/(h^2))*(theta(i+1)+2*theta(i)-theta(i-1))+(B/(2*h))*(v(i+1)-v(i-1)); % internal nodes for 0<z<d

end

% v equation

% for i = N: corresponds to i=9-5 =4 =N

rhsode(N+2,1) = C/(h^2)*(v(N-1) + 2*v(N-2) - v(N-3))+ ...

+ G/(h^3)*(-theta(N-3) +3*theta(N-2) - 3*theta(N-1) + theta(N))...

+xi/(2*h)*(theta(N-1) - theta(N-3)); %involves special third derivative

% for i = 8, this corresponds i = 8-5 =3

% let j =N+4-(N+1) = 3

for j=3:N-1

theta(N+1) = 0;

v(N+1) = 0;

rhsode(j,1)= C/(h^2)*(v(j+1) + 2*v(j) - v(j-1))+ ...

+ G/(2*h^3)*(v(j+2) -2*theta(j+1) + 2*theta(j-1) - theta(j-2))...

+ xi/(2*h)*(theta(j+1) - theta(j-1)); %uses central difference for third derivative for theta

end

% for i = N: corresponds to i=9-5 =4 =N

theta(N+1) = 0;

v(N+1) = 0;

rhsode(2*N,1) = C/(h^2)*(v(N+1) - 2*v(N) + 2*v(N-1))+ ...

+ G/(h^3)*(-theta(N-2) +3*theta(N-1) -3*theta(N) + theta(N+1))...

+xi/(2*h)*(theta(N+1) - theta(N-1)); %involves special third derivative

end