function?richards
%?Solution?of?the?Richards?equation?????
%????using?MATLAB?pdepe???????????????????
%
%???$Ekkehard?Holzbecher??$Date:?2006/07/13?$
%
%?Soil?data?for?Guelph?Loam?（Hornberger?and?Wiberg?2005）
%?m-file?based?partially?on?a?first?version?by?Janek?Greskowiak?
%--------------------------------------------------------------------------
L?=?60;???????????????????%?length?[cm]
T?=?4;????????????????????%?maximum?time?[h]
qr?=?0.102;???????????????%?residual?water?content?殘余含水量
f?=?0.368;????????????????%?porosity??飽和含水量
s1=0.4;???????????????????%下滲率
s2=0;?????????????????????%底部壓力水頭
a?=?0.0335;???????????????%?van?Genuchten?parameter?[1/L]
n?=?2;????????????????????%?van?Genuchten?parameter
m?=?1-1/n;
ks?=0.00922;????????????????%?saturated?conductivity?[L/T]
N=100;
M=25;
x?=?linspace（0LN）;
t?=?linspace（0TM）;?

options=odeset（‘RelTol‘1e-4‘AbsTol‘1e-4‘NormControl‘‘off‘‘InitialStep‘1e-7）
u?=?pdepe（0@unsatpde@unsatic@unsatbcxtoptionss1s2qrfanks）;


figure;???????????%表示壓力水頭（基質勢）分別在t=0h4h24h時刻隨土壤深度分布圖
size（u）
plot?（xu（1:）‘-‘）;
xlabel（‘Depth?[L]‘）;
ylabel（‘Pressure?Head?[L]‘）;
legend（‘t=0h‘‘t=4h‘‘t=24h‘）;
title（‘壓力水頭分別在t=0h4h24h時刻隨土壤深度分布圖‘）;

figure;??????????????????????%表示壓力水頭分別在z=0cm16cm36cm60cm位置處隨時間變化分布圖
plot?（tu（:1）‘-‘）;
xlabel（‘T?[h]‘）;
ylabel（‘Pressure?Head?[L]‘）;
legend（‘z=0cm‘‘z=16cm‘‘z=36cm‘‘z=60cm‘）
title（‘壓力水頭分別在z=0cm16cm36cm60cm位置處隨時間變化分布圖‘）;


for?j=1:length（t）
????for?i=1:length（x）
?