Solve a complex differential equation over different intervals

Hello,

I am trying to solve this third order differential equation whose constants have different values depending on the intervals considered (3). I wrote a piece of code and tried some things especially with "switch case" but I don't have the necessary knowledge to solve correctly so I gave up, if someone could help me it would be very kind.

Here is the differential equation and the code

clear all
clc

etilde=6/25;

zi=-1;
zf=1;
phii=1;
phif=-1;

n=10^4;

yinit = [phii phif];
init = bvpinit([linspace(zi,zi+etilde,n) linspace(zi+etilde,zf-etilde,n) linspace(zf-etilde,1,n)],yinit);
sol = bvp4c(@(x,y)f(x,y), @bc, init);
x=[linspace(zi,zi+etilde,n) linspace(zi+etilde,zf-etilde,n) linspace(zf-etilde,1,n)];
y = deval(sol, x);
plot(x,y(1,:));
legend('phi(z)')
title('Potentiel électrique en fonction de l épaisseur adimensionnée')
xlabel({'z'})
ylabel('phi(z)')

function dphidz = f(x,y) % equations being solved

epsilon=10^-5;

A1=[51.8 19.8 51.8];
A2=[0 3.35 0];
A3=epsilon*[0.499 1.58 0.499];
A4=[1.18E-1 1.09E-2 1.18E-1];
A5=[2.44 0.617 2.44];
A6=0.6;
Phi2m=[1 0.83 1];
etilde=6/25;

y = zeros(3,1);

dphidz = zeros(3,1);

dphidz(1)=y(2);
dphidz(2)=y(3);

if (-1<z)&&(z<-1+etilde)
     % z in [-1 -1+etilde]
    j=1;
    dphidz(3)=A6/(A3(j)*A1(j)*Phi2m(j))*y(2)*(A1(j)-A2(j)+A3(j)*y(3))/...
    (A4(j)+A5(j)/A6*(1-A1(j)*Phi2m(j)/(A1(j)-A2(j)+A3(j)*y(3))));
elseif (-1+etilde<z)&&(z<1-etilde)   % z in [-1+etilde 1-etilde]
    j=2;
	dphidz(3)=A6/(A3(j)*A1(j)*Phi2m(j))*y(2)*(A1(j)-A2(j)+A3(j)*y(3))/...
	(A4(j)+A5(j)/A6*(1-A1(j)*Phi2m(j)/(A1(j)-A2(j)+A3(j)*y(3))));
else    % z in [1-etilde 1]
	j=3;
	dphidz(3)=A6/(A3(j)*A1(j)*Phi2m(j))*y(2)*(A1(j)-A2(j)+A3(j)*y(3))/...
	(A4(j)+A5(j)/A6*(1-A1(j)*Phi2m(j)/(A1(j)-A2(j)+A3(j)*y(3))));
end

end

function res = bc(phiL,phiR)

epsilon=10^-5;

A3=epsilon*[0.499 1.58 0.499];

res = [phiL(1,1) - 1              
       phiR(1,1) - phiL(1,2)    
       A3(1)*phiR(2,1) - A3(2)*phiL(2,2)
       A3(1)*phiL(2,1) - A3(2)*phiL(2,3)
       A3(1)*phiR(3,1) - A3(2)*phiL(3,2)
       phiR(1,2) - phiL(1,3)     
       A3(2)*phiR(2,2) - A3(3)*phiL(2,3)
       A3(2)*phiR(3,2) - A3(3)*phiL(3,3)
       phiR(1,3) + 1];       
end

Answer 1

[darova]

Here are some notes:

• there is no need of creating special mesh. You have 3d order diff equation so you need 3 yinit

yinit = [phii phif 1];
% init = bvpinit([linspace(zi,zi+etilde,n) linspace(zi+etilde,zf-etilde,n) linspace(zf-etilde,1,n)],yinit);
init = bvpinit(-1:.1:1,yinit);

• This part of the code can be shorter and more readable

• too many boundary contions. You need only 3

% res = [phiL(1,1) - 1              
%        phiR(1,1) - phiL(1,2)    
%        A3(1)*phiR(2,1) - A3(2)*phiL(2,2)
%        A3(1)*phiL(2,1) - A3(2)*phiL(2,3)
%        A3(1)*phiR(3,1) - A3(2)*phiL(3,2)
%        phiR(1,2) - phiL(1,3)     
%        A3(2)*phiR(2,2) - A3(3)*phiL(2,3)
%        A3(2)*phiR(3,2) - A3(3)*phiL(3,3)
%        phiR(1,3) + 1];       
res = [phiL(1) - 1
    phiR(1) + 1
    phiR(2) - phiL(2)];