Sum_{j=1}^{4} D_{1j} * (drond^{2}z_{j} / drond x^{2} + drond^{2}z_{j} / drond y^{2}) + R_{1} = 0

Sum_{j=1}^{4} D_{2j} * (drond^{2}z_{j} / drond x^{2} + drond^{2}z_{j} / drond y^{2}) + R_{2} = 0

Sum_{j=1}^{4} D_{3j} * (drond^{2}z_{j} / drond x^{2} + drond^{2}z_{j} / drond y^{2}) + R_{3} = 0

Sum_{j=1}^{4} D_{4j} * (drond^{2}z_{j} / drond x^{2} + drond^{2}z_{j} / drond y^{2}) + R_{4} = 0

Boundary conditions: x = 0,1 z_{1} = 1 ; x = 0,1 z_{2,3,4} = 0

So, I have a system with four partial diffrential equations. and I want to plot the solution for every equations.

What it's wrong in my code ? The plots of solutions must be a surface (almost like a bell), but is not.

Code:

```
clc;
clear;
close all;
a=0;b=1;n=5;
h = (b-a)/(n+1); % step
I=eye((n-1)*(n-1));
A11=zeros((n-1)*(n-1));A12=zeros((n-1)*(n-1));A13=zeros((n-1)*(n-1)); A14=zeros((n-1)*(n-1));
A21=zeros((n-1)*(n-1));A22=zeros((n-1)*(n-1));A23=zeros((n-1)*(n-1));A24=zeros((n-1)*(n-1));
A31=zeros((n-1)*(n-1));A32=zeros((n-1)*(n-1));A33=zeros((n-1)*(n-1));A34=zeros((n-1)*(n-1));
A41=zeros((n-1)*(n-1));A42=zeros((n-1)*(n-1));A43=zeros((n-1)*(n-1));A44=zeros((n-1)*(n-1));
b1=zeros((n-1)*(n-1),1);b2=zeros((n-1)*(n-1),1);b3=zeros((n-1)*(n-1),1);b4=zeros((n-1)*(n-1),1);
k1=1;k2=0.4;k3=0.5;k4=0.2;k5=0.1;
% diffusion matrix
D11=1; D12=-0.43; D13=-0.19; D14=-0.265;
D21=-0.14; D22=1.475; D23=-0.025; D24=0.235;
D31=0.105; D32=0.07; D33=1.475; D34=-0.19;
D41=-0.085; D42=-0.245; D43=0.27; D44=1.17;
[x,y] = meshgrid(a:h:b); % uniform mesh
% Matrix A
K1D = spdiags(ones(n,1)*[-1 2 -1],-1:1,n-1,n-1);
I1D = speye(size(K1D));
Delta = kron(K1D,I1D)+kron(I1D,K1D);
% A1i
A11=D11*Delta-(k1+k3)*I*h^2;
A12=D12*Delta;
A13=D13*Delta;
A14=D14*Delta-k5*I*h^2;
%A2i
A21=D21*Delta-k1*I*h^2;
A22=D22*Delta+k2*I*h^2;
A23=D23*Delta-k4*I*h^2;
A24=D24*Delta;
%A3i
A31=D31*Delta;
A32=D32*Delta-k2*I*h^2;
A33=D33*Delta+k4*I*h^2;
A34=D34*Delta;
%A4i
A41=D41*Delta-k3*I*h^2;
A42=D42*Delta;
A43=D43*Delta;
A44=D44*Delta+k5*I*h^2;
% b
b1(1)=(-1/h^2)*D11;
b2(1)=(-1/h^2)*D21;
b3(1)=(-1/h^2)*D31;
b4(1)=(-1/h^2)*D41;
b1(n-1)=(-1/h^2)*D11;
b2(n-1)=(-1/h^2)*D21;
b3(n-1)=(-1/h^2)*D31;
b4(n-1)=(-1/h^2)*D41;
% aprox solution
A=[A11 A12 A13 A14;A21 A22 A23 A24;A31 A32 A33 A34;A41 A42 A43 A44];
bn=[b1;b2;b3;b4];
ub=reshape(bn,(n-1)*(n-1)*4,1);
v=inv(A)*ub;
% split the solution vector
pas=size(A11);
v1=v(1:n+2);
v2=v(pas+1:2*pas);
v3=v(2*pas+1:3*pas);
v4=v(3*pas+1:4*pas);
vv1=zeros(n+2,n+2);
vv2=zeros(n+2,n+2);
vv3=zeros(n+2,n+2);
vv4=zeros(n+2,n+2);
% boundary conditions
vv1=[ones(n+2,1) v1 ones(n+2,1) ones(n+2,1) ones(n+2,1) ones(n+2,1) ones(n+2,1)];
vv1=[zeros(n+2,1) v1 zeros(n+2,1) zeros(n+2,1) zeros(n+2,1) zeros(n+2,1) zeros(n+2,1)];
vv1=[zeros(n+2,1) v1 zeros(n+2,1) zeros(n+2,1) zeros(n+2,1) zeros(n+2,1) zeros(n+2,1)];
vv1=[zeros(n+2,1) v1 zeros(n+2,1) zeros(n+2,1) zeros(n+2,1) zeros(n+2,1) zeros(n+2,1)];
% plot
figure
set(gcf,'DefaultAxesFontSize',8,'PaperPosition', [0 0 3.5 3.5])
subplot(2,2,1);
surf(x,y,vv1);
title('Numerical solution for Ecuation1');
subplot(2,2,2);
surf(x,y,vv2);
title('Numerical solution for Ecuation2');
subplot(2,2,3);
surf(x,y,vv3)
title('Numerical solution for Ecuation3');
subplot(2,2,4);
surf(x,y,vv4);
title('Numerical solution for Ecuation4');
```