Matlab程式程式碼：

clear s t x0 y0 a b
syms s t
f=sin(s)*sin(t);
a=16.0;
b=16.0;
x0=-8;
y0=-8;
m=40;
n=40;
h1=a/m;
h2=b/n;
hold  % 將新舊圖形保持在同一個圖中
for i=1:m   % i 從1每次以1為距離差一直加到m
    s=x0+(i-1)*h1;  
    for j=1:n  % j 從1每次以1為距離差一直加到n
        t=y0+(j-1)*h2;
        d=eval(f);  % 執行式子f
        y1=t+2/3*h1*d;
        if abs(y1-t)>2/3*h2  % 檢查(y1-t)的絕對值是否大於2/3*h2
            x1=s+1/d*h2*2/3;  % 如果上面的if語句為真則執行
            plot([s,x1],[t,t+h2*2/3]) % 以向量"[s,x1]"為x軸，以向量“[t,t+h2*2/3]”為y軸繪圖
        else
            plot([s,s+h1*2/3],[t,y1]) % 以向量"[s,s+h1*2/3]"為x軸，以向量“[t,y1]”為y軸繪圖 
        end
    end
end
title('dy/dx=sinx*siny');  % 以“dy/dx=sinx*siny”為該圖形的標題
xlabel('x');  % 標記x軸
ylabel('y')   % 標記y軸
 

執行結果：