How to solve system of PDE's using crank nicolson method to get graphical interpretations of equations? How to get skin friction and nusselt number using this code?

i have attached my matlb code for solving system of PDE ∂u/∂t=(∂^2 u)/〖∂η〗^2 +α (∂^3 u)/〖∂η〗^3 +δcφ+Grθ-Mu

∂θ/∂t=1/Pr (∂^2 θ)/〖∂η〗^2 +D_f (∂^2 φ)/〖∂η〗^2

∂φ/∂t=1/S_c (∂^2 φ)/〖∂η〗^2 +S_r (∂^2 θ)/〖∂η〗^2

when η→0

when η→0:u=sin(wt),θ=1 ,φ=1, u^'=1

when η→∞:u→0,θ→0,φ→0

But,i'm confused how to set boundary conditions in this code. Also i want to set these skin friction and nusselt number expressions: C_f=-(∂u/∂η) at η=0

Nu〖Re〗_x^(-1)=-(∂θ/∂η) at η=0

in this code;

function yy= crank

n=105;

h=0.01;

m=0.1;

B=0; %B=0.0001;

Gr=2; %Gr=0.1=Delt. Gs=Dels

A=2; %A=alpha

Pr=5;

Gs=1; %3.5

Sc=1.5;

M=9; %M=30,12

Df=5;

Sr=10;

Ec=0.1; % 0.01-0.1(0-1)

t(1)=0;

s(1)=0;

for j=2:n

t(j)=t(j-1)+h;

end

for j=2:n

s(j)=s(j-1)+m;

end

k=1;

for j=1:n

aa(j,k)=exp(-0.15*Pr^2*t(j)*Ec^-0.1/Gr^0.1); % Guess M*Gr*Gs*S

bb(j,k)=exp(-t(j)*Sc/(Ec)^0.01*M); %bb(j,k)=exp(-Pr^2*t(j)/Ec);

cc(j,k)=exp(-Sc*2*t(j)/Ec^0.3);

end

a{1,k}=[1 0 0;0 1 0;0 0 1];

for j=2:n-1

a{j,k}=[-(1/h)-(1/h^2)-2*A-M Gr Gs;0 -(1/h)-(1/Pr*(h^2)) -Df/h^2;0 -Sr/h^2 -(1/h)-(1/Sc*h^2)];

end

a{n,k}=[1 0 0;0 1 0;0 0 1];

for j=2:n-1

b{j,k}=[(1/2*h^2)+(A/h^3) 0 0;0 1/2*Pr*h^2 Df/2*h^2;0 Sr/2*h^2 1/2*Sc*h^2]; % Lower diagonal entries

end

b{n,k}=[0 0 0;0 0 0;0 0 0];

c{1,k}=[0 0 0;0 0 0;0 0 0];

for j=2:n-1

c{j,k}=[(1/2*h^2)+(A/h^3) 0 0;0 1/2*Pr*h^2 Df/2*h^2;0 Sr/2*h^2 1/2*Sc*h^2];

end

r1(1,k)=0;

r2(1,k)=0;

r3(1,k)=0;

for j=2:n-1

r1(j,k)=-(exp(t(j))/h)-(1/2*h^2)*(1)*(exp(-t(j+1))-2*exp(-t(j))+exp(-t(j-1)))-Gr*exp(-t(j))-Gs*exp(t(j))+M*exp(-t(j))-(A/2*h^3)*(-exp(-t(j+1))+2*exp(-t(j))-exp(-t(j-1)))...

+(aa(j,k)/h)-(1/2*h^2)*(1)*(aa(j+1,k)-2*aa(j,k)+aa(j-1,k))-(A/2*h^3)*( aa(j+1,k)^2-2*aa(j+1,k)*aa(j,k))-Gr*bb(j,k)-Gs*cc(j,k)+M*aa(j,k);

r2(j,k)=(exp((-t(j)))/h)-(1/2*Pr*h^2)*(exp(-t(j+1))-2* exp(-t(j))+exp(-t(j-1)))-Df*(1/2*h^2)*(exp(-t(j+1))-2* exp(-t(j))...

+exp(-t(j-1)))+(bb(j,k)/h)-(1/2*Pr*h^2)*(bb(j+1,k)-2*bb(j,k)+bb(j-1,k))-(Df/2*h^2)*(cc(j+1,k)-2*cc(j,k)+cc(j-1,k));

r3(j,k)=-(exp((-t(j)))/h)-(1/2*Sc*h^2)*(exp(-t(j+1))-2*exp(-t(j))+exp(-t(j-1)))-(Sr/2*h^2)*(exp(-t(j+1))-2*exp(-t(j))+exp(-t(j-1)))+ (cc(j,k)/h)...

-(1/2*Sc*h^2)*(cc(j+1,k)-2*cc(j,k)+cc(j-1,k))-(Sr/2*h^2)*(bb(j+1,k)-2*bb(j,k)+bb(j-1,k));

end

r1(n,k)= 0;

r2(n,k)=0;

r3(n,k)=0;

for j=1:n

rr{j,k}=[r1(j,k);r2(j,k);r3(j,k)];

end

gamma{1,k}=inv(a{1,k})*c{1,k};

for j=2:n-1

a{j,k}=a{j,k}-(b{j,k}*gamma{j-1,k});

gamma{j,k}=inv(a{j,k})*c{j,k};

end

y{1}=inv(a{1})*rr{1};

for j=2:n

y{j}=inv(a{j})*(rr{j}-b{j}*y{j-1});

end

x{n}=y{n};

for j=n-1:-1:1

x{j}=y{j}-(gamma{j})*x{j+1};

end

for j=n:-1:1

u(j,k)=x{j}(1,1);

z(j,k)=x{j}(2,1);

q(j,k)=x{j}(3,1);

end

for j=1:n

xx(j,k)= aa(j,k)+u(j,k);

yy(j,k)= bb(j,k)+z(j,k);

zz(j,k)= cc(j,k)+q(j,k);

end

for j=1:n-1

%U(j)=-((1+(1/B))*((xx(j+1)-xx(j))/h))

%U(j)=-((yy(j+1)-yy(j))/2*h)

end

%ee=meshgrid(xx);

%V=trapz(ee);

%eee=meshgrid(V);

%[XX,YY]=meshgrid(t,s);

%XXX=meshgrid(xx);

%surf(XX,YY,XXX)

%contour(eee)

figure(1)

plot(t,xx,'LineWidth',2,'MarkerSize',16,'linestyle','-','color','k')

xlabel('\eta','Interpreter','tex','FontSize',16,'FontWeight','bold');

ylabel('u','Interpreter','tex','FontSize',16,'FontWeight','bold');

grid on

hold on

figure(2)

plot(t,yy,'LineWidth',2,'MarkerSize',16,'linestyle','-','color','k')

xlabel('\eta','Interpreter','tex','FontSize',16,'FontWeight','bold');

ylabel('\Theta(\eta)','Interpreter','tex','FontSize',16,'FontWeight','bold');

grid on

hold on

figure(3)

plot(t,zz,'LineWidth',2,'MarkerSize',16,'linestyle','--','color','r')

xlabel('\eta','Interpreter','tex','FontSize',16,'FontWeight','bold');

ylabel('\Phi(\eta)','Interpreter','tex','FontSize',16,'FontWeight','bold');

grid on

hold on

[Torsten]

I might be mistaken, but the third derivative of u with respect to eta suggests that you need 7 boundary conditions in total, not only 6.

[Sharqa]

yes i mistakenly forget one boundary condition. Now, I added missing one.can you help me how to fix these boundary condition in given code