I wrote a following program. Everything is ok to my knowledge. But it shows the error; "Index exceeds matrix dimensions."
clear all
clc
format longEng
syms x y z kh
a=[4;0.5;70];
% The Newton-Raphson iterations starts here
KH=linspace(0,0.55,12);
h=4;
q=100;
gma=14; nq=2*q/(gma*h);
delta=10;
phi=30;
lam=0.1;
A=lam*nq/(1+nq);
kv=0;
% kh=0;
da1=delta*(pi/180); da2=-delta*(pi/180); pha1=phi*(pi/180); pha2=phi*(pi/180);
dp1=delta*(pi/180); dp2=delta*(pi/180); php1=phi*(pi/180); php2=phi*(pi/180);
psi=atan(kh/(1-kv));
a=pha1+da1;
b=pha1-psi;
c=psi+da1;
alphac=atan((sin(a)*sin(b)+(sin(a)^2+sin(b)^2+sin(a)*cos(a)*sin(b)*cos(b)+A*cos(c)*cos(a)*sin(b))^0.5)/(A*cos(c)+sin(a)*cos(b)));
w=0.5*gma*h^2*(1/tan(alphac));
B=h*((1/tan(alphac))-lam);
ka1=((1+nq)*(1-A*tan(alphac))*(cos(b)-(sin(b)/tan(alphac))))/(cos(psi)*(cos(a)+tan(alphac)*sin(a)));
va2=asin(sin(da2)/sin(pha2))-asin(sin(psi)/sin(pha2))-da2-psi;
ka2=(1/cos(psi))*(cos(da2)*((cos(da2)-sqrt(sin(pha2)^2-sin(da2)^2)))/(cos(psi)+sqrt(sin(pha2)^2-sin(psi)^2)))*exp(-va2*tan(pha2));
vp1=asin(sin(dp1)/sin(php1))+asin(-sin(psi)/sin(php1))+dp1+psi;
kp1=(1/cos(psi))*(cos(dp1)*((cos(dp1)+sqrt(sin(php1)^2-sin(dp1)^2)))/(cos(psi)-sqrt(sin(php1)^2-sin(psi)^2)))*exp(vp1*tan(php1));
vp2=asin(sin(dp2)/sin(php2))+asin(-sin(psi)/sin(php2))+dp2+psi;
kp2=(1/cos(psi))*(cos(dp2)*((cos(dp2)+sqrt(sin(php2)^2-sin(dp2)^2)))/(cos(psi)-sqrt(sin(php2)^2-sin(psi)^2)))*exp(vp2*tan(php2));
% Qa1=q*ka1*(h+x); Qp2=q*kp2*y;
sinda1=sin(da1); sindp1=sin(dp1); sinda2=-sin(da2); sindp2=sin(dp2);
cosda1=cos(da1); cosdp1=cos(dp1); cosda2=cos(da2); cosdp2=cos(dp2);
% pa1=ka1*gma*0.5*(h+x)^2;
pp1=kp1*gma*0.5*(x^2);
pa1=0.5*gma*(1-kv)*ka1*h^2;
pa2=ka2*gma*(x*y+0.5*(y^2)); pp2=kp2*gma*(y*(h+x)+(0.5*(y^2)));
zp1=x/3;
za1=((((w*(1-kv)-pa1*sin(da1))*cos(pha1)*h)/(sin(alphac)*cos(alphac-pha1)))-(q*B*(lam*h+0.5*B)))/(pa1*cos(da1));
zp2=((0.5*(h+x)*(y^2))+((y^3)/3))/(((h+x)*y)+(0.5*(y^2)));
za2=((0.5*x*(y^2))+((y^3)/3))/((x*y)+(0.5*(y^2)));
e2=(pp1*cosdp1)+(pa2*cosda2)-(pa1*cosda1)-(pp2*cosdp2);
e3=(pp1*cosdp1*zp1)+(pp2*cosdp2*zp2)-(pa1*cosda1*za1)-(pa2*cosda2*za2);
g=[e2; e3]
J=jacobian([e2, e3], [x, y]);
A=zeros(2,numel(KH));
for i=1:numel(KH)
del=1;
indx=0;
while del>1e-6
gnum = vpa(subs(g,[x,y,kh],[a(1),a(2),KH(i)]));
Jnum = vpa(subs(J,[x,y,kh],[a(1),a(2),KH(i)]));
delx = -Jnum\gnum;
a = a + delx;
del = max(abs(gnum));
indx = indx + 1;
end
A(:,i)=double(a)
end
