求教各位高手帮忙指点程序

gaussmathematic

小弟用有限差分法隐式格式编了一个程序求解偏微分方程，但在调试时一直有问题，运行非常缓慢，出来结果也不可用。调试了很长时间都没弄好，无奈之下想到论坛，所以把程序贴上来，想请各位帮忙看看。万分感谢！
程序如下：

function Ure=myfun(T,G,g) %由于两个方向上步长一致时，计算非线性方程组是对称的，只能用T2来区别到底是哪个方向采用隐式或显示
%for i=2:20
%UU(i)=T(i,:);
%end
f1=0.2;
f2=0.01;
UU=T;%[4.2 4.20 4.2 4.2 4.2 4.2 4.2 4.2 4.2 4.2 4.2 4.2 4.2 4.2 4.2 4.2 4.2 4.2 4.2 4.2 4.2];
gg=g;
h=1e-4;%0.01;
k=1e-9;%0.0001;
r=k/(h*h);
r1=k/h;
%G=14e+6;
GG=G*k;
%x1=4.2;  %此条件限定了每层边界都为4.2K 采用偏导数条件
x11=4.2  %此条件限定了每层边界都为4.2K
theta=3045/12800;
%C=f1*6550*(0.152+2.10*10^(-3)*T^3)+f2*6380*(1.1-0.4*T+5e-2*T^2-4e-5*T^3)+(1-f1-f2)*8940*(7.582e-4*T^3);
%K=f1*(0.38887*T^0.1532-0.371)+f2*(-3.5332+9.6273*T-0.1282*T^2+5e-4*T^3)+(1-f1-f2)*(10^(1.9153+1.7202*log(T)-0.4146*log(T)^2));
%******syms符号变量，要用20个字符
%for i=2:20
%ff(i)=r*A(U(i)+U(i-1))*(U(i)-U(i-1))+B(U(i))*(U(i)-UU(i))-r*A(U(i)+U(i+1))*(U(i+1)-U(i));
%end
xx1=0.0001;
for i=1:10
    x(i)=(i-1)*h;
end
%************为了利用matlab中diff函数计算倒数****************
syms  x1 x2 x3 x4 x5 x6 x7 x8 x9 x10
C1=f1*6550*(0.152+2.10*10^(-3)*x1^3)+f2*6380*(1.1-0.4*x1+5e-2*x1^2-4e-5*x1^3)+(1-f1-f2)*8940*(7.582e-4*x1^3);
C2=f1*6550*(0.152+2.10*10^(-3)*x2^3)+f2*6380*(1.1-0.4*x2+5e-2*x2^2-4e-5*x2^3)+(1-f1-f2)*8940*(7.582e-4*x2^3);
C3=f1*6550*(0.152+2.10*10^(-3)*x3^3)+f2*6380*(1.1-0.4*x3+5e-2*x3^2-4e-5*x3^3)+(1-f1-f2)*8940*(7.582e-4*x3^3);
C4=f1*6550*(0.152+2.10*10^(-3)*x4^3)+f2*6380*(1.1-0.4*x4+5e-2*x4^2-4e-5*x4^3)+(1-f1-f2)*8940*(7.582e-4*x4^3);
C5=f1*6550*(0.152+2.10*10^(-3)*x5^3)+f2*6380*(1.1-0.4*x5+5e-2*x5^2-4e-5*x5^3)+(1-f1-f2)*8940*(7.582e-4*x5^3);
C6=f1*6550*(0.152+2.10*10^(-3)*x6^3)+f2*6380*(1.1-0.4*x6+5e-2*x6^2-4e-5*x6^3)+(1-f1-f2)*8940*(7.582e-4*x6^3);
C7=f1*6550*(0.152+2.10*10^(-3)*x7^3)+f2*6380*(1.1-0.4*x7+5e-2*x7^2-4e-5*x7^3)+(1-f1-f2)*8940*(7.582e-4*x7^3);
C8=f1*6550*(0.152+2.10*10^(-3)*x8^3)+f2*6380*(1.1-0.4*x8+5e-2*x8^2-4e-5*x8^3)+(1-f1-f2)*8940*(7.582e-4*x8^3);
C9=f1*6550*(0.152+2.10*10^(-3)*x9^3)+f2*6380*(1.1-0.4*x9+5e-2*x9^2-4e-5*x9^3)+(1-f1-f2)*8940*(7.582e-4*x9^3);
C10=f1*6550*(0.152+2.10*10^(-3)*x10^3)+f2*6380*(1.1-0.4*x10+5e-2*x10^2-4e-5*x10^3)+(1-f1-f2)*8940*(7.582e-4*x10^3);
%C11=f1*6550*(0.152+2.10*10^(-3)*x11^3)+f2*6380*(1.1-0.4*x11+5e-2*x11^2-4e-5*x11^3)+(1-f1-f2)*8940*(7.582e-4*x11^3);
%C12=f1*6550*(0.152+2.10*10^(-3)*x12^3)+f2*6380*(1.1-0.4*x12+5e-2*x12^2-4e-5*x12^3)+(1-f1-f2)*8940*(7.582e-4*x12^3);
%C13=f1*6550*(0.152+2.10*10^(-3)*x13^3)+f2*6380*(1.1-0.4*x13+5e-2*x13^2-4e-5*x13^3)+(1-f1-f2)*8940*(7.582e-4*x13^3);
%C14=f1*6550*(0.152+2.10*10^(-3)*x14^3)+f2*6380*(1.1-0.4*x14+5e-2*x14^2-4e-5*x14^3)+(1-f1-f2)*8940*(7.582e-4*x14^3);
%C15=f1*6550*(0.152+2.10*10^(-3)*x15^3)+f2*6380*(1.1-0.4*x15+5e-2*x15^2-4e-5*x15^3)+(1-f1-f2)*8940*(7.582e-4*x15^3);
%C16=f1*6550*(0.152+2.10*10^(-3)*x16^3)+f2*6380*(1.1-0.4*x16+5e-2*x16^2-4e-5*x16^3)+(1-f1-f2)*8940*(7.582e-4*x16^3);
%C17=f1*6550*(0.152+2.10*10^(-3)*x17^3)+f2*6380*(1.1-0.4*x17+5e-2*x17^2-4e-5*x17^3)+(1-f1-f2)*8940*(7.582e-4*x17^3);
%C18=f1*6550*(0.152+2.10*10^(-3)*x18^3)+f2*6380*(1.1-0.4*x18+5e-2*x18^2-4e-5*x18^3)+(1-f1-f2)*8940*(7.582e-4*x18^3);
%C19=f1*6550*(0.152+2.10*10^(-3)*x19^3)+f2*6380*(1.1-0.4*x19+5e-2*x19^2-4e-5*x19^3)+(1-f1-f2)*8940*(7.582e-4*x19^3);
%C20=f1*6550*(0.152+2.10*10^(-3)*x20^3)+f2*6380*(1.1-0.4*x20+5e-2*x20^2-4e-5*x20^3)+(1-f1-f2)*8940*(7.582e-4*x20^3);

K1=f1*(0.38887*(x2/2+x1/2)^0.1532-0.371)+f2*(-3.5332+9.6273*(x2/2+x1/2)-0.1282*(x2/2+x1/2)^2+5e-4*(x2/2+x1/2)^3)+(1-f1-f2)*(10^(1.9153+1.7202*log10((x2/2+x1/2))-0.4146*log10((x2/2+x1/2))^2));
K2=f1*(0.38887*(x3/2+x2/2)^0.1532-0.371)+f2*(-3.5332+9.6273*(x3/2+x2/2)-0.1282*(x3/2+x2/2)^2+5e-4*(x3/2+x2/2)^3)+(1-f1-f2)*(10^(1.9153+1.7202*log10((x3/2+x2/2))-0.4146*log10((x3/2+x2/2))^2));
K3=f1*(0.38887*(x4/2+x3/2)^0.1532-0.371)+f2*(-3.5332+9.6273*(x4/2+x3/2)-0.1282*(x4/2+x3/2)^2+5e-4*(x4/2+x3/2)^3)+(1-f1-f2)*(10^(1.9153+1.7202*log10((x4/2+x3/2))-0.4146*log10((x4/2+x3/2))^2));
K4=f1*(0.38887*(x5/2+x4/2)^0.1532-0.371)+f2*(-3.5332+9.6273*(x5/2+x4/2)-0.1282*(x5/2+x4/2)^2+5e-4*(x5/2+x4/2)^3)+(1-f1-f2)*(10^(1.9153+1.7202*log10((x5/2+x4/2))-0.4146*log10((x5/2+x4/2))^2));
K5=f1*(0.38887*(x6/2+x5/2)^0.1532-0.371)+f2*(-3.5332+9.6273*(x6/2+x5/2)-0.1282*(x6/2+x5/2)^2+5e-4*(x6/2+x5/2)^3)+(1-f1-f2)*(10^(1.9153+1.7202*log10((x6/2+x5/2))-0.4146*log10((x6/2+x5/2))^2));
K6=f1*(0.38887*(x7/2+x6/2)^0.1532-0.371)+f2*(-3.5332+9.6273*(x7/2+x6/2)-0.1282*(x7/2+x6/2)^2+5e-4*(x7/2+x6/2)^3)+(1-f1-f2)*(10^(1.9153+1.7202*log10((x7/2+x6/2))-0.4146*log10((x7/2+x6/2))^2));
K7=f1*(0.38887*(x8/2+x7/2)^0.1532-0.371)+f2*(-3.5332+9.6273*(x8/2+x7/2)-0.1282*(x8/2+x7/2)^2+5e-4*(x8/2+x7/2)^3)+(1-f1-f2)*(10^(1.9153+1.7202*log10((x8/2+x7/2))-0.4146*log10((x8/2+x7/2))^2));
K8=f1*(0.38887*(x9/2+x8/2)^0.1532-0.371)+f2*(-3.5332+9.6273*(x9/2+x8/2)-0.1282*(x9/2+x8/2)^2+5e-4*(x9/2+x8/2)^3)+(1-f1-f2)*(10^(1.9153+1.7202*log10((x9/2+x8/2))-0.4146*log10((x9/2+x8/2))^2));
K9=f1*(0.38887*(x10/2+x9/2)^0.1532-0.371)+f2*(-3.5332+9.6273*(x10/2+x9/2)-0.1282*(x10/2+x9/2)^2+5e-4*(x10/2+x9/2)^3)+(1-f1-f2)*(10^(1.9153+1.7202*log10((x10/2+x9/2))-0.4146*log10((x10/2+x9/2))^2));
K10=f1*(0.38887*(x11/2+x10/2)^0.1532-0.371)+f2*(-3.5332+9.6273*(x11/2+x10/2)-0.1282*(x11/2+x10/2)^2+5e-4*(x11/2+x10/2)^3)+(1-f1-f2)*(10^(1.9153+1.7202*log10((x11/2+x10/2))-0.4146*log10((x11/2+x10/2))^2));
%K11=f1*(0.38887*(x12/2+x11/2)^0.1532-0.371)+f2*(-3.5332+9.6273*(x12/2+x11/2)-0.1282*(x12/2+x11/2)^2+5e-4*(x12/2+x11/2)^3)+(1-f1-f2)*(10^(1.9153+1.7202*log10((x12/2+x11/2))-0.4146*log10((x12/2+x11/2))^2));
%K12=f1*(0.38887*(x13/2+x12/2)^0.1532-0.371)+f2*(-3.5332+9.6273*(x13/2+x12/2)-0.1282*(x13/2+x12/2)^2+5e-4*(x13/2+x12/2)^3)+(1-f1-f2)*(10^(1.9153+1.7202*log10((x13/2+x12/2))-0.4146*log10((x13/2+x12/2))^2));
%K13=f1*(0.38887*(x14/2+x13/2)^0.1532-0.371)+f2*(-3.5332+9.6273*(x14/2+x13/2)-0.1282*(x14/2+x13/2)^2+5e-4*(x14/2+x13/2)^3)+(1-f1-f2)*(10^(1.9153+1.7202*log10((x14/2+x13/2))-0.4146*log10((x14/2+x13/2))^2));
%K14=f1*(0.38887*(x15/2+x14/2)^0.1532-0.371)+f2*(-3.5332+9.6273*(x15/2+x14/2)-0.1282*(x15/2+x14/2)^2+5e-4*(x15/2+x14/2)^3)+(1-f1-f2)*(10^(1.9153+1.7202*log10((x15/2+x14/2))-0.4146*log10((x15/2+x14/2))^2));
%K15=f1*(0.38887*(x16/2+x15/2)^0.1532-0.371)+f2*(-3.5332+9.6273*(x16/2+x15/2)-0.1282*(x16/2+x15/2)^2+5e-4*(x16/2+x15/2)^3)+(1-f1-f2)*(10^(1.9153+1.7202*log10((x16/2+x15/2))-0.4146*log10((x16/2+x15/2))^2));
%K16=f1*(0.38887*(x17/2+x16/2)^0.1532-0.371)+f2*(-3.5332+9.6273*(x17/2+x16/2)-0.1282*(x17/2+x16/2)^2+5e-4*(x17/2+x16/2)^3)+(1-f1-f2)*(10^(1.9153+1.7202*log10((x17/2+x16/2))-0.4146*log10((x17/2+x16/2))^2));
%K17=f1*(0.38887*(x18/2+x17/2)^0.1532-0.371)+f2*(-3.5332+9.6273*(x18/2+x17/2)-0.1282*(x18/2+x17/2)^2+5e-4*(x18/2+x17/2)^3)+(1-f1-f2)*(10^(1.9153+1.7202*log10((x18/2+x17/2))-0.4146*log10((x18/2+x17/2))^2));
%K18=f1*(0.38887*(x19/2+x18/2)^0.1532-0.371)+f2*(-3.5332+9.6273*(x19/2+x18/2)-0.1282*(x19/2+x18/2)^2+5e-4*(x19/2+x18/2)^3)+(1-f1-f2)*(10^(1.9153+1.7202*log10((x19/2+x18/2))-0.4146*log10((x19/2+x18/2))^2));
%K19=f1*(0.38887*(x20/2+x19/2)^0.1532-0.371)+f2*(-3.5332+9.6273*(x20/2+x19/2)-0.1282*(x20/2+x19/2)^2+5e-4*(x20/2+x19/2)^3)+(1-f1-f2)*(10^(1.9153+1.7202*log10((x20/2+x19/2))-0.4146*log10((x20/2+x19/2))^2));
%K20=f1*(0.38887*(x21/2+x20/2)^0.1532-0.371)+f2*(-3.5332+9.6273*(x21/2+x20/2)-0.1282*(x21/2+x20/2)^2+5e-4*(x21/2+x20/2)^3)+(1-f1-f2)*(10^(1.9153+1.7202*log10((x21/2+x20/2))-0.4146*log10((x21/2+x20/2))^2));




gre1=subs(-2*r*K1*(x2-x1)/(theta-1),{x1,x2},{gg(1,1),gg(1,2)});
gre2=subs(r*K1*(x2-x1)/(theta-1)-r*K2*(x3-x2)/(theta-1),{x1,x2,x3},{gg(2,1),gg(2,2),gg(2,3)});
gre3=subs(r*K2*(x3-x2)/(theta-1)-r*K3*(x4-x3)/(theta-1),{x2,x3,x4},{gg(3,1),gg(3,2),gg(3,3)});
gre4=subs(r*K3*(x4-x3)/(theta-1)-r*K4*(x5-x4)/(theta-1),{x3,x4,x5},{gg(4,1),gg(4,2),gg(4,3)});
gre5=subs(r*K4*(x5-x4)/(theta-1)-r*K5*(x6-x5)/(theta-1),{x4,x5,x6},{gg(5,1),gg(5,2),gg(5,3)});
gre6=subs(r*K5*(x6-x5)/(theta-1)-r*K6*(x7-x6)/(theta-1),{x5,x6,x7},{gg(6,1),gg(6,2),gg(6,3)});
gre7=subs(r*K6*(x7-x6)/(theta-1)-r*K7*(x8-x7)/(theta-1),{x6,x7,x8},{gg(7,1),gg(7,2),gg(7,3)});
gre8=subs(r*K7*(x8-x7)/(theta-1)-r*K8*(x9-x8)/(theta-1),{x7,x8,x9},{gg(8,1),gg(8,2),gg(8,3)});
gre9=subs(r*K8*(x9-x8)/(theta-1)-r*K9*(x10-x9)/(theta-1),{x8,x9,x10},{gg(9,1),gg(9,2),gg(9,3)});
gre10=subs(r*K9*(x10-x9)/(theta-1)-r*K10*(x11-x10)/(theta-1),{x9,x10,x11},{gg(10,1),gg(10,2),gg(10,3)});
%gre11=subs(r*K10*(x11-x10)-r*K11*(x12-x11),{x10,x11,x12},{gg(11,1),gg(11,2),gg(11,3)});
%gre12=subs(r*K11*(x12-x11)-r*K12*(x13-x12),{x11,x12,x13},{gg(12,1),gg(12,2),gg(12,3)});
%gre13=subs(r*K12*(x13-x12)-r*K13*(x14-x13),{x12,x13,x14},{gg(13,1),gg(13,2),gg(13,3)});
%gre14=subs(r*K13*(x14-x13)-r*K14*(x15-x14),{x13,x14,x15},{gg(14,1),gg(14,2),gg(14,3)});
%gre15=subs(r*K14*(x15-x14)-r*K15*(x16-x15),{x14,x15,x16},{gg(15,1),gg(15,2),gg(15,3)});
%gre16=subs(r*K15*(x16-x15)-r*K16*(x17-x16),{x15,x16,x17},{gg(16,1),gg(16,2),gg(16,3)});
%gre17=subs(r*K16*(x17-x16)-r*K17*(x18-x17),{x16,x17,x18},{gg(17,1),gg(17,2),gg(17,3)});
%gre18=subs(r*K17*(x18-x17)-r*K18*(x19-x18),{x17,x18,x19},{gg(18,1),gg(18,2),gg(18,3)});
%gre19=subs(r*K18*(x19-x18)-r*K19*(x20-x19),{x18,x19,x20},{gg(19,1),gg(19,2),gg(19,3)});
%gre20=subs(r*K19*(x20-x19)-r*K20*(x21-x20),{x19,x20},{gg(20,1),gg(20,2)});

ff1=-2*r*K1*(x2-x1)+C1*(x1-UU(1))-1/xx1*r1*K1*(x2-x1)-GG+gre1;%注意格式前正负号
ff2=r*K1*(x2-x1)+C2*(x2-UU(2))-r*K2*(x3-x2)-1/x(2)*r1*K2*(x3-x2)-GG+gre2;
ff3=r*K2*(x3-x2)+C3*(x3-UU(3))-r*K3*(x4-x3)-1/x(3)*r1*K3*(x4-x3)-GG+gre3;
ff4=r*K3*(x4-x3)+C4*(x4-UU(4))-r*K4*(x5-x4)-1/x(4)*r1*K4*(x5-x4)-GG+gre4;
ff5=r*K4*(x5-x4)+C5*(x5-UU(5))-r*K5*(x6-x5)-1/x(5)*r1*K5*(x6-x5)-GG+gre5;
ff6=r*K5*(x6-x5)+C6*(x6-UU(6))-r*K6*(x7-x6)-1/x(6)*r1*K6*(x7-x6)-GG+gre6;
ff7=r*K6*(x7-x6)+C7*(x7-UU(7))-r*K7*(x8-x7)-1/x(7)*r1*K7*(x8-x7)-GG+gre7;
ff8=r*K7*(x8-x7)+C8*(x8-UU(8))-r*K8*(x9-x8)-1/x(8)*r1*K8*(x9-x8)-GG+gre8;
ff9=r*K8*(x9-x8)+C9*(x9-UU(9))-r*K9*(x10-x9)-1/x(9)*r1*K9*(x10-x9)-GG+gre9;
ff10=r*K9*(x10-x9)+C10*(x10-UU(10))-r*K10*(x11-x10)-1/x(10)*r1*K10*(x11-x10)-GG+gre10;
%ff11=r*K10*(x11-x10)+C11*(x11-UU(11))-r*K11*(x12-x11)-1/x(11)*r1*K11*(x12-x11)+gre11;%-GG;
%ff12=r*K11*(x12-x11)+C12*(x12-UU(12))-r*K12*(x13-x12)-1/x(12)*r1*K12*(x13-x12)+gre12;%-GG;
%ff13=r*K12*(x13-x12)+C13*(x13-UU(13))-r*K13*(x14-x13)-1/x(13)*r1*K13*(x14-x13)+gre13;%-GG;
%ff14=r*K13*(x14-x13)+C14*(x14-UU(14))-r*K14*(x15-x14)-1/x(14)*r1*K14*(x15-x14)+gre14;%-GG;
%ff15=r*K14*(x15-x14)+C15*(x15-UU(15))-r*K15*(x16-x15)-1/x(15)*r1*K15*(x16-x15)+gre15;%-GG;
%ff16=r*K15*(x16-x15)+C16*(x16-UU(16))-r*K16*(x17-x16)-1/x(16)*r1*K16*(x17-x16)+gre16;%-GG;
%ff17=r*K16*(x17-x16)+C17*(x17-UU(17))-r*K17*(x18-x17)-1/x(17)*r1*K17*(x18-x17)+gre17;%-GG;
%ff18=r*K17*(x18-x17)+C18*(x18-UU(18))-r*K18*(x19-x18)-1/x(18)*r1*K18*(x19-x18)+gre18;%-GG;
%ff19=r*K18*(x19-x18)+C19*(x19-UU(19))-r*K19*(x20-x19)-1/x(19)*r1*K19*(x20-x19)+gre19;%-GG;
%ff20=r*K19*(x20-x19)+C20*(x20-UU(20))-r*K20*(x21-x20)-1/x(20)*r1*K20*(x21-x20)+gre20;%-GG;
ff=[ff1 ff2 ff3 ff4 ff5 ff6 ff7 ff8 ff9 ff10];% ff11 ff12 ff13 ff14 ff15 ff16 ff17 ff18 ff19 ff20]';

%G=[diff(ff(1),x1) diff(ff(1),x2) 0 0
%   diff(ff(2),x1) diff(ff(2),x2) diff(ff(2),x3) 0
%   0 diff(ff(3),x2) diff(ff(3),x3) diff(ff(3),x4)
%   0 0 diff(ff(4),x3) diff(ff(4),x4)];

%
%   
G=[diff(ff(1),x1) diff(ff(1),x2) 0 0 0 0 0 0 0 0
    diff(ff(2),x1) diff(ff(2),x2) diff(ff(2),x3) 0 0 0 0 0 0 0
    0 diff(ff(3),x2) diff(ff(3),x3) diff(ff(3),x4) 0 0 0 0 0 0
    0 0 diff(ff(4),x3) diff(ff(4),x4) diff(ff(4),x5) 0 0 0 0 0
    0 0 0 diff(ff(5),x4) diff(ff(5),x5) diff(ff(5),x6) 0 0 0 0
    0 0 0 0 diff(ff(6),x5) diff(ff(6),x6) diff(ff(6),x7) 0 0 0
    0 0 0 0 0 diff(ff(7),x6) diff(ff(7),x7) diff(ff(7),x8) 0 0
    0 0 0 0 0 0 diff(ff(8),x7) diff(ff(8),x8) diff(ff(8),x9) 0
    0 0 0 0 0 0 0 diff(ff(9),x8) diff(ff(9),x9) diff(ff(9),x10)
    0 0 0 0 0 0 0 0 diff(ff(10),x9) diff(ff(10),x10)];
   U0=[15.18 15.19 15.20 15.21 15.22 15.23 15.24 15.25 15.26 15.27]';
   ff0=subs(ff,{x1 x2 x3 x4 x5 x6 x7 x8 x9 x10},{U0(1) U0(2) U0(3) U0(4) U0(5) U0(6) U0(7) U0(8) U0(9) U0(10)})';
   G0=subs(G,{x1 x2 x3 x4 x5 x6 x7 x8 x9 x10},{U0(1) U0(2) U0(3) U0(4) U0(5) U0(6) U0(7) U0(8) U0(9) U0(10)});  
   
  
% G=[diff(ff(1),x1) diff(ff(1),x2) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
%   diff(ff(2),x1) diff(ff(2),x2) diff(ff(2),x3) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
%    0 diff(ff(3),x2) diff(ff(3),x3) diff(ff(3),x4) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
%    0 0 diff(ff(4),x3) diff(ff(4),x4) diff(ff(4),x5) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
%    0 0 0 diff(ff(5),x4) diff(ff(5),x5) diff(ff(5),x6) 0 0 0 0 0 0 0 0 0 0 0 0 0 0
  %  0 0 0 0 diff(ff(6),x5) diff(ff(6),x6) diff(ff(6),x7) 0 0 0 0 0 0 0 0 0 0 0 0 0
%    0 0 0 0 0 diff(ff(7),x6) diff(ff(7),x7) diff(ff(7),x8) 0 0 0 0 0 0 0 0 0 0 0 0
%   0 0 0 0 0 0 diff(ff(8),x7) diff(ff(8),x8) diff(ff(8),x9) 0 0 0 0 0 0 0 0 0 0 0
%    0 0 0 0 0 0 0 diff(ff(9),x8) diff(ff(9),x9) diff(ff(9),x10) 0 0 0 0 0 0 0 0 0 0
  %  0 0 0 0 0 0 0 0 diff(ff(10),x9) diff(ff(10),x10) diff(ff(10),x11) 0 0 0 0 0 0 0 0 0
%    0 0 0 0 0 0 0 0 0 diff(ff(11),x10) diff(ff(11),x11) diff(ff(11),x12) 0 0 0 0 0 0 0 0
  %  0 0 0 0 0 0 0 0 0 0 diff(ff(12),x11) diff(ff(12),x12) diff(ff(12),x13) 0 0 0 0 0 0 0
%   0 0 0 0 0 0 0 0 0 0 0 diff(ff(13),x12) diff(ff(13),x13) diff(ff(13),x14) 0 0 0 0 0 0
%    0 0 0 0 0 0 0 0 0 0 0 0 diff(ff(14),x13) diff(ff(14),x14) diff(ff(14),x15) 0 0 0 0 0
%    0 0 0 0 0 0 0 0 0 0 0 0 0 diff(ff(15),x14) diff(ff(15),x15) diff(ff(15),x16) 0 0 0 0
%   0 0 0 0 0 0 0 0 0 0 0 0 0 0 diff(ff(16),x15) diff(ff(16),x16) diff(ff(16),x17) 0 0 0
%   0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diff(ff(17),x16) diff(ff(17),x17) diff(ff(17),x18) 0 0
%   0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diff(ff(18),x17) diff(ff(18),x18) diff(ff(18),x19) 0
%   0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diff(ff(19),x18) diff(ff(19),x19) diff(ff(19),x20)
%   0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diff(ff(20),x19) diff(ff(20),x20)];
%ff';
% G;
%***************把符号表达式转换成变量表达式用来赋初值U0***************
%U0=[4.19 4.20 4.21 4.22 4.23 4.24 4.25 4.26 4.27 4.28 4.29 4.30 4.31 4.32 4.33 4.34 4.35 4.36 4.37]';
%U0=[15.18 15.19 15.20 15.21 15.22 15.23 15.24 15.25 15.26 15.27 15.28 15.29 15.30 15.31 15.32 15.33 15.34 15.35 15.36 15.37]';
%U0=10*[15.18 15.19 15.20 15.21 15.22 15.23 15.24 15.25 15.26 15.27 15.28 15.29 15.30 15.31 15.32 15.33 15.34 15.35 15.36 15.37]';
%ff0=subs(ff,{x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 x17 x18 x19 x20},{U0(1) U0(2) U0(3) U0(4) U0(5) U0(6) U0(7) U0(8) U0(9) U0(10) U0(11) U0(12) U0(13) U0(14) U0(15) U0(16) U0(17) U0(18) U0(19) U0(20)})';
%G0=subs(G,{x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 x17 x18 x19 x20},{U0(1) U0(2) U0(3) U0(4) U0(5) U0(6) U0(7) U0(8) U0(9) U0(10) U0(11) U0(12) U0(13) U0(14) U0(15) U0(16) U0(17) U0(18) U0(19) U0(20)});



%************ 用fsolve求解**********



%***********牛顿法求解**********
U1=U0-inv(G0)*ff0;
while(norm(U1-U0)>=1.0e-15)
   U0=U1;
   ff0=subs(ff,{x1 x2 x3 x4 x5 x6 x7 x8 x9 x10},{U0(1) U0(2) U0(3) U0(4) U0(5) U0(6) U0(7) U0(8) U0(9) U0(10)})';
   G0=subs(G,{x1 x2 x3 x4 x5 x6 x7 x8 x9 x10},{U0(1) U0(2) U0(3) U0(4) U0(5) U0(6) U0(7) U0(8) U0(9) U0(10)});  
    U1=U0-inv(G0)*ff0;
end



U1;
Ure=[U1;4.2];



主程序差分格式求解
function myfun
clear all
f1=0.2 ;     
f2=0.01;

G=14e+6;%+8.6300e+6;
h=1e-4;%0.01;
k=1e-9;
r1=k/(h*h);
r2=k/h;

%A=K/c
%G=14e+6;
%Q=8.6300e+6;
T=[];
x=[];
y=[];  
TT=[];

%const=[GG GG GG GG GG GG 0 0 0 0 0 0 0 0 0 0 0 0 0 0]';
dP=[];%修改0426
g=[];
g1=[];
flag=0;
flag1=0;
theta=3045/12800;
%初始条件
for i=1:11%41     %i=1;i<=20;i++                     %不能用m必须代入实际的数值
   for j=1:31              %时间第一层，在时间为0时，即起始
x(i)=(i-1)*h; %空间步长
y(i)=(i-1)*h;
T(i,j,1)=4.2;
end
end

for j=1:21
    z(j)=(j-1)*k;
end

for j=1:21                     %指定边界条件；两方向上都认为恒温
    for l=1:11
        T(l,31,j)=4.2;
    end
end

for j=1:21                   %指定边界条件；两方向上都认为恒温
    for m=1:31
        T(11,m,j)=4.2;
    end
end
%边界条件
%for j=1:101
   %T(1,j)=4.2;              %T对x的导数边界条件为0
%   T(22,j)=4.2;               %分为20等分，但总层数为21，包含起始层
                               %空间第一层，即表面；
%   y(j)=(j-1)*k;
%end

%C=f1*6550*(0.152+2.10*10^(-3)*T^3)+f2*6380*(1.1-0.4*T+5e-2*T^2-4e-5*T^3)+(1-f1-f2)*8940*(7.582e-4*T^3);
%K=f1*(0.38887*T^0.1532-0.371)+f2*(-3.5332+9.6273*T-0.1282*T^2+5e-4*T^3)+(1-f1-f2)*(10^(1.9153+1.7202*log(T)-0.4146*log(T)^2));
%******************用牛顿法求解非线性格式******************

for j=1:5
    G1=G;
    clear maplemex
     for l=1:10
            g(l,:)=[T(l,2,j),T(l,1,j),T(l,2,j)];%[4.2 4.2 4.2];        % 计算l向隐式时，需要m向的三个值
        end
T(:,1,j+1)=twodijizuobiaonewton10(T(:,1,j),G1,g);
  
    for m=2:30
        flag=flag+1
          %if flag==8
          %      for l=1:4
          %  g(l,:)=[(T(l,m-1,j)+T(l,m,j))/2,(T(l,m,j)+T(l,m+1,j))/2,(T(l,m+1,j)+T(l,m-1,j))/2]
         %       end
         % else
        for l=1:10
            g(l,:)=[theta*T(l,m,j)+(1-theta)*T(l,m+1,j),theta*T(l,m+1,j)+(1-theta)*T(l,m-1,j),theta*T(l,m-1,j)+(1-theta)*T(l,m,j)];%[T(l,m-1,j),T(l,m,j),T(l,m+1,j)]   %[4.2 4.2 4.2];     % 计算l向隐式时，需要m向的三个值
        end
      
        T(:,m,j+1)=twodijizuobiaonewton10(T(:,m,j),G1,g); %flag中gre2和gre3值不匹配
       %   end
    end
    g;
    shuchu=T(:,3,3);
  end  


  
  
  
  

for j=6:20
    G2=0;
    clear maplemex
       for l=1:10
            g(l,:)=[T(l,2,j),T(l,1,j),T(l,2,j)]; %[4.2 4.2 4.2];      % 计算l向隐式时，需要m向的三个值
        end
  T(:,1,j+1)=twodijizuobiaonewton10(T(:,1,j),G2,g);
  
   
   for m=2:30
       flag1=flag1+1
        for l=1:10
            g(l,:)=[theta*T(l,m,j)+(1-theta)*T(l,m+1,j),theta*T(l,m+1,j)+(1-theta)*T(l,m-1,j),theta*T(l,m-1,j)+(1-theta)*T(l,m,j)];%[T(l,m-1,j),T(l,m,j),T(l,m+1,j)];%[4.2 4.2 4.2];
        end
        T(:,m,j+1)=twodijizuobiaonewton10(T(:,m,j),G2,g)';
    end

end   
   




T
for i=1:21
    Tp2(i)=T(2,2,i);
    Tp4(i)=T(4,4,i);
end
plot(z,Tp2,'go')
hold on
plot(z,Tp4,'ro')
% hold on
% plot(z,T(6,6,:),'g*')
% hold on
% plot(z,T(8,8,:),'r*')
% hold on
% plot(z,T(10,10,:),'g+')

gaussmathematic

自己先顶一下，大侠们，施施援手啊。快崩溃了。

07047224

你太牛啦，这么多代码，看着都晕，不会，顶你下

mwglikai

可否把你的程序所应用的背景贴出来。。

gaussmathematic

mwglikai 发表于 2012-6-25 17:00
可否把你的程序所应用的背景贴出来。。

几何模型比较简单，就是一个圆柱体，圆柱体由三层构成，由内到外依次为第一种非线性材料，金属铜，第二种非线性材料。这里的非线性是指这两种材料的电阻是非线性的。也就是说整个导体是三种材料复合而成。
现在要考察给整个圆柱导体加直流电流或交流电流后，分析导体的温度分布，总体上是一个热电耦合场分析，但是其中电阻非线性，想用ansys求解，问了一下好像ansys不能做非线性电阻的电场分析，所以采用matlab编了程序求解，但是运行得不到结果且容易出现矩阵奇异这样的错误。

liuyang5299

对矩阵G0求逆时要注意G0是否可逆