通过subs求解这两千多个函数的过程，速度实在太慢了

如下边这段代码，是matlab通过subs求解两千多个函数的过程，速度实在太慢了，请看看，有什么好的办法可以提高计算速度？最好附上代码，多谢！

clear all;clc;
tic
dx=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

没有好的想法，希望其他大佬帮下忙。

你这边将其分解成离散的点进行所以相对而言比较慢，数据多的话建议采用拟合，然后利用相关的算法求近似解