valavanje sferaste membrane -- > simulacija

Pozdravljeni,
mogoče bi mi vi lahko pomagali. Namreč sem se lotil simulirat valovanje na sferični membrani. Seveda prvo sem se lotil simulirati ravne membrane, kjer tudi deluje super s formulo:

U(i,j,k)=r*(U(i-1,j,k-1)+U(i+1,j,k-1)+U(i,j-1,k-1)+U(i,j+1,k-1))-U(i,j,k-2)+(2-4*r)*U(i,j,k-1);

pri čemer je i in j je točka na membrani, k pa je naslednjo stanje, r je hitrost širjenje valovanja, pa mora biti manjše kot 0.5.

recimo kos programa je tako če bo hotel kdo od vas programirat:

for k=3:M
U(end/2-2:end/2+2,2,2)=0.1*sin(2*pi*k/30);
for i=1:m
for j=1:n
if (D(i,j)==0)
U(i,j,3)=r*(U(i-1,j,2)+U(i+1,j,2)+U(i,j-1,2)+U(i,j+1,2))-U(i,j,1)+(2-4*r)*U(i,j,2);
end
end
end
end

OK vse lepo in prav, ampak jaz sem se lotil da bi naredil valovanje za sferično membrano (valovanje vode po sferi) pa sem prišel do take formule

U(i,j,k)=r/cos(theta(j))^2*(U(i-1,j,k-1)+U(i+1,j,k-1)+cos(theta(j))^2*(U(i,j-1,k-1)+U(i,j+1,k-1)))-U(i,j,k-2)+(2-2*r*(1/cos(theta(j))^2+1))*U(i,j,k-1);

i=phi, j=theta
delta phi = delta theta, tako da je prišla zgornja formula malo bolj poenostavljena in to pomeni da imamo 2 več i-jev kot j-jev

saj na začetku celo deluje ampak čez nekaj časa pa vse podvija (divergira). A bi znal kdo kaj pomagat, zakaj tako podivja?

Poišči kje je napaka prevelika. Potem se pa lahko igraš s stvarmi kot je razvijanje vsega kar leze in gre v Taylorjevo vrsto (predvsem v upanju, da se kaj pokrajša in/ali izpostavi). Če tako ne gre, je, resda grda, rešitev, da strateško postaviš kak pogoj, ki vlovi očitno napačne rezultate. Pri tem si lahko pomagaš tudi s tem, da hraniš nekaj zadnjih rezultatov (ali pa kar vse in potem zaradi pridobljene hitrosti izvajanja uporabiš bolj fino delitev).