syms z11(x) z12(x) z21(x) z22(x) z31(x) z32(x)
eqn1 = diff(z11) == z12;
eqn2 = diff(z12) + z12 + z22 + z31 - 5 == 0;
eqn3 = diff(z21) == z22;
eqn4 = diff(z22) + z12 + z22 + z32 - 4 == 0;
eqn5 = diff(z31) == z32;
eqn6 = diff(z32) - z22 + z11 + z31 - 3 == 0;
eqns = [eqn1,eqn2,eqn3,eqn4,eqn5,eqn6];
conds = [z11(0)==0,z21(0)==0,z31(0)==0,z11(1)==0,z21(1)==0,z31(1)==0];
sol = dsolve(eqns,conds)
bvpfcn = @(x,y)[y(2);-y(2)-y(4)-y(5)+5;y(4);-y(2)-y(4)-y(6)+4;y(6);y(4)-y(1)-y(5)+3];
bcfcn = @(ya,yb) [ya(1);ya(3);ya(5);yb(1);yb(3);yb(5)];
xmesh = linspace(0,1,100);
solinit = bvpinit(xmesh, [0 0 0 0 0 0]);
sol = bvp4c(bvpfcn, bcfcn, solinit);
plot(sol.x, [sol.y(1,:);sol.y(3,:);sol.y(5,:)])
