solution by newton raphson method for nonlinear equations

I am writing a code for solving two non linear simultaneous equations using newton raphson method. I am not able to link the g and J for different variables with newton raphson method. As I am new to matlab. Please help and thank in advance.

dbstop if error
clear all
clc
format longEng
syms x y h
phi=(pi/180)*39;
delta=(2*phi)/3;
gma=18.4;
a=[1.5;0.632];
% kh=0.3;
H=linspace(0.5,4,8);
% h=4;
lam=0;
q=50;
nq=2*q/(gma*(h+x));
A=lam*nq/(1+nq);
kh=0;
kv=0;
psi=atan(kh/(1-kv));
beta=1.2;
alfa=1.2;
R1=-1; 
R3=-1;
delm1=0.5*(1-R1)*delta;
delm3=-0.5*(1-R3)*delta;
m=phi+delm1;
b=phi-psi;
c=psi+delm1;

alphac=atan((sin(m)*sin(b)+(sin(m)^2*sin(b)^2+sin(m)*cos(m)*sin(b)*cos(b)+A*cos(c)*cos(m)*sin(b))^0.5)/(A*cos(c)+sin(m)*cos(b)));

kg=(tan(alphac-phi)+(kh/(1-kv)))/(tan(alphac)*(cos(delm1)+sin(delm1)*tan(alphac-phi)));
r=1-lam*tan(alphac);
kq=r*kg;

pg=0.5*gma*(1-kv)*kg*(h+x)^2*cos(delm1);
pq=(1-kv)*q*kq*(h+x);

k3=(2*cos(phi-psi)^2)/(cos(phi-psi)^2*(1+R3)+cos(psi)*cos(delm3+psi)*(1-R3)*(1+sqrt((sin(phi+delm3)*sin(phi-psi))/cos(delm3+psi)))^2);


y1=0:0.01:1;
y2=0:0.01:1;
mmm=size(y1);
mm=mmm(1,2);

for i=1:mm
    if (y1(i)>=1-(1/beta)&& y1(i)<=1)
        R2(i)=3*(beta*(1-y1(i))).^0.5;
    else
        R2(i)=3;
    end
    if (R2(i)>=0 && R2(i)<1)
        delm2(i)=0.5*(1-R2(i)).*delta;
        k2(i)=(2*cos(phi-psi)^2)/(cos(phi-psi)^2*(1+R2(i))+cos(psi)*cos(delm2(i)+psi)...
            *(1-R2(i))*(1+sqrt((sin(phi+delm2(i))*sin(phi-psi))/cos(delm2(i)+psi)))^2);
    else 
        delm2(i)=0.5*(R2(i)-1)*delta;
        k2(i)=1+0.5*(R2(i)-1).*((cos(phi-psi).^2./(cos(psi).*(cos(delm2(i)+psi).*...
            (-sqrt((sin(phi+delm2(i)).*sin(phi-psi))./(cos(delm2(i)+psi)))+1).^2)))-1);
    end
    if (y2(i)>=0 && y2(i)<=(1/alfa))
        R4(i)=3*(alfa*(y2(i)))^0.5;
    else
        R4(i)=3;
    end
     if (R4(i)>=0 && R4(i)<=1)
         delm4(i)=0.5*(1-R4(i)).*delta;
         k4(i)=(2*cos(phi-psi)^2)/(cos(phi-psi)^2*(1+R4(i))+cos(psi)*cos(delm4(i)+psi)...
             *(1-R4(i))*(1+sqrt((sin(phi+delm4(i))*sin(phi-psi))/cos(delm4(i)+psi)))^2);
     else
         delm4(i)=0.5*(R4(i)-1)*delta;
         k4(i)=1+0.5*(R4(i)-1)*((cos(phi-psi)^2/(cos(psi)*(cos(delm4(i)+psi)*(-sqrt((sin(phi+delm4(i))...
             *sin(phi-psi))/(cos(delm4(i)+psi)))+1)^2)))-1);
     end
     
h2(i)=gma.*x.^2.*k2(i).*y1(i).*cos(delm2(i));
h4(i)=gma.*y.^2.*k4(i).*cos(delm4(i)).*y2(i)+gma.*(x+h).*y.*k4(i).*cos(delm4(i));
H3_a(i)=k3.*y2(i).*cos(delm3);
% H3_b=matlabFunction(k3*cos(delm3));
h3=gma.*y.^2.*H3_a+gma.*x.*y.*k3.*cos(delm3);%integral(H3_b,0,1);
HF=h2-h4+h3-pg-pq;
%
M2(i)=k2(i).*y1(i).*cos(delm2(i)).*(1-y1(i));
m2=gma.*x.^3.*M2;
M4_a(i)=y2(i).^2;
M4_b(i)=y2(i);
m4=gma.*y.^3.*k4(i).*cos(delm4(i)).*M4_a+gma.*(x+h).*y.^2.*k4(i).*cos(delm4(i)).*M4_b;
M3_a(i)=k3.*y2(i).^2.*cos(delm3);
M3_b(i)=k3.*y2(i).*cos(delm3);
m3=gma.*y.^3.*M3_a+gma.*x.*y.^2.*M3_b;
MF=m2+m4-m3-pg.*(h+x).*(1/3)-0.5.*pq.*(h+x);
g=[HF; MF];
J=jacobian([HF, MF], [x, y]);
end

% Z=zeros(2,numel(H));
for j=1:numel(H)
    del=1;
    indx=0;

while del>1e-6
   gnum = vpa(subs(g,[x,y],[a(1),a(2)]));
   Jnum = vpa(subs(J,[x,y],[a(1),a(2)]));
   delx = -Jnum.\gnum;
   a = a + delx;
   del = max(abs(gnum));
   indx = indx + 1;
end
 Z(:,j)=double(a)
end

'NEWTON-RAPHSON SOLUTION CONVERGES IN ITERATIONS',indx,
'FINAL VALUES OF a ARE';a,

[Walter Roberson]

http://www.mathworks.co.uk/matlabcentral/answers/43073-a-guide-to-tags

[Walter Roberson]

What is the difference between this and the earlier https://www.mathworks.com/matlabcentral/answers/470367-error-using-symengine-the-dimensions-of-matrices-or-vectors-are-incompatible ?