The process of solving these more than 2000 functions through subs is too slow

As shown in the following code , yes matlab adopt subs The process of solving more than 2000 functions , The speed is too slow , Please have a look , Is there any good way to improve the calculation speed ？ It's better to attach the code , Thank you very much ！

clear all;clc;ticdx=0.1;   dy=0.1;   dt=0.1;nx=100;    ny=100;a=1;      b=1;P=sym('p_',[nx ny]);xx1=-(diag(ones(1,nx),0)-diag(ones(1,nx-1),1)).*dx;xx2=(diag(ones(1,nx),0)-diag(ones(1,nx-1),-1)).*dx; yy1=-(diag(ones(1,ny),0)-diag(ones(1,ny-1),-1)).*dy;yy2=(diag(ones(1,ny),0)-diag(ones(1,ny-1),1)).*dy; XX1=xx1*P;      XX2=xx2*P;YY1=P*yy1;      YY2=P*yy2; fx1=(1-a*exp(-b*abs(XX1./dx))).*(XX1./dx);fx2=(1-a*exp(-b*abs(XX2./dx))).*(XX2./dx);fy1=(1-a*exp(-b*abs(YY1./dy))).*(YY1./dx);fy2=(1-a*exp(-b*abs(YY2./dy))).*(YY2./dx); fx=(fx1-fx2)./dx;fy=(fy1-fy2)./dy; fxy=fx+fy;A=ones(nx,ny);ft=0.5.*(P-A)/dt;F=fxy-ft;F1=reshape(F,1,nx*ny);P1=reshape(P,1,nx*ny);r=ones(1,nx*ny);Fxyt=double(subs(F1,P1,r));toc