[ÇóÖú] ÓÃMATLABÇó½â΢·Ö·½³Ì×éʱPDEPEº¯ÊýÖв»ÖªµÀÔõô±íʾ·½³Ì×é

function [c,f,s] = pdex3pde(x,t,u,DuDx)
%±»µ÷Óõĺ¯Êý£¬Æä×÷ÓÃÊÇÃèÊöƫ΢·Ö·½³Ì
c = [1; 1; 1; 0];
f =   
s =   

function pdex44 %pdex3ΪÎļþÃû
clc
m = 0;
x = [0 0.005 0.01 0.05];
t = [0 0.005 0.01 0.05];

sol = pdepe(m,@pdex44pde,@pdex44ic,@pdex44bc,x,t);
%pdepeº¯Êý£¬ÓÃÓÚÖ±½ÓÇó½âƫ΢·Ö·½³Ì£¬ÆäÐÎʽΪsol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)
%ÏÂÃæΪ×÷ͼ²½Öè

u1 = sol(:,:,1);
% ui = sol(:,:,i) is an approximation to the ith component of the solution vector u .
u2 = sol(:,:,2);

figure
surf(x,t,u1)
title('u1(x,t)')
xlabel('Distance x')
ylabel('Time t')

figure
surf(x,t,u2)
title('u2(x,t)')
xlabel('Distance x')
ylabel('Time t')

% --------------------------------------------------------------
function [c,f,s] = pdex44pde(x,t,u,DuDx)
%±»µ÷Óõĺ¯Êý£¬Æä×÷ÓÃÊÇÃèÊöƫ΢·Ö·½³Ì
c = [1; 1; 1; 0];
f = [80*DuDx(1); (2.4*10^6)*DuDx(2)-6*(10^5)*u(2).*DuDx(4); (2.4*10^6)*DuDx(3)-(6*10^5)*u(3).*DuDx(4); DuDx(4)];
s = [2888*u(1).*DuDx(1).*DuDx(4)-(5.2*10^-5)*u(1).*(u(1)-u(2))+10350*u(2); 10350*u(2); (9.6*10^-14)*u(2).*((DuDx(4)).^2)-(3.8*10^-13)*DuDx(2).*DuDx(4)-(2.099*10^-9)*u(3); (1.8*10^-8)*(u(1)-u(2))];

% --------------------------------------------------------------
function u0 = pdex44ic(x);
%´Ëº¯ÊýÊÇÓÃÀ´ÃèÊö³õÖµ
u0 = [10^6; 10^7; 0; 0];
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex44bc(xl,ul,xr,ur,t)
%´Ëº¯ÊýÓÃÀ´ÃèÊö±ß½çÌõ¼þ
pl = [0; -(1772036.5/4)*ul(2); (1772036.5/4)*ul(3); -500];
ql = [1; 1; 1; 1];
pr = [0; -(1772036.5/4)*ul(2); (1772036.5/4)*ul(3); 0];
qr = [1; 1; 1; 1];