Well, my computer is a little slower than yours, but this is definitely better! Now (with Mat_proc_newer), the creation of the matrix is (on my computer) approximatelly as fast as in MS Excel and the minimization process lasts around 26s which is significantly faster (than in Excel).

Accuracy: the results are the same (as in Excel), so it's OK.

Well, my computer is a little slower than yours, but this is definitely better! Now (with Mat_proc_newer), the creation of the matrix is (on my computer) approximatelly as fast as in MS Excel and the minimization process lasts around 26s which is significantly faster (than in Excel).

Accuracy: the results are the same (as in Excel), so it's OK.

@reemeaaaah You assign A[1]:=B[3] for k=3, but then you call for every k this:

sol[k] := dsolve({ic[k], subs(A[1] = B[1], A[2] = B[2], ode)}, U(t), numeric)

which "throws away" the assignment you made before.

You can type e.g.:

for k to 4 do         
  tk := 365*k;     
  V := rhs(sol[k-1](tk)[2]);     
  ic[k] := U(tk) = sigma*V;
    if k = 3 then sol[k] := dsolve({ic[k], subs(A[1] = B[3], A[2] = B[2], ode)}, U(t), numeric)
                 else sol[k] := dsolve({ic[k], subs(A[1] = B[1], A[2] = B[2], ode)}, U(t), numeric)  
    end if;   
end do;

@reemeaaaah You assign A[1]:=B[3] for k=3, but then you call for every k this:

sol[k] := dsolve({ic[k], subs(A[1] = B[1], A[2] = B[2], ode)}, U(t), numeric)

which "throws away" the assignment you made before.

You can type e.g.:

for k to 4 do         
  tk := 365*k;     
  V := rhs(sol[k-1](tk)[2]);     
  ic[k] := U(tk) = sigma*V;
    if k = 3 then sol[k] := dsolve({ic[k], subs(A[1] = B[3], A[2] = B[2], ode)}, U(t), numeric)
                 else sol[k] := dsolve({ic[k], subs(A[1] = B[1], A[2] = B[2], ode)}, U(t), numeric)  
    end if;   
end do;

1) You also don't need this:

sol := proc(rkf45_t) end;

2) There are two square brackets missing in the plots[odeplot] command, so you should type e.g.:

plots[odeplot](sol, [[t, T1(t)], [t, T2(t), color = blue]], 0 .. 20);

1) You also don't need this:

sol := proc(rkf45_t) end;

2) There are two square brackets missing in the plots[odeplot] command, so you should type e.g.:

plots[odeplot](sol, [[t, T1(t)], [t, T2(t), color = blue]], 0 .. 20);

@Red Horse I am sorry I don't know if I understood correctly your question.

A(n) is the "function" representing the area of the n-th fractal (I used the word square, but from n=2 on it is not a square, of course).

A(1)=1, since the area of the square of side-length equal to 1 is 1^2=1. Area of the n-th fractal is the area of the previous one + the areas of the newly arisen squares. The number of these new squares is equal to N(n-1), since there is one new square at every side of the previous fractal. The area of one such new square is L(n)^2, since the length of one side of this square is L(n).

So this

A(1)=1:
A(n)=A(n-1)+N(n-1)*L(n)^2:

can be the answer. But this is a recurrence which can be solved, so we can get A(n) as a function of n which is more practical. That's why I used the ?rsolve command just to get a different (more practical) formula for the same thing.

@Red Horse I am sorry I don't know if I understood correctly your question.

A(n) is the "function" representing the area of the n-th fractal (I used the word square, but from n=2 on it is not a square, of course).

A(1)=1, since the area of the square of side-length equal to 1 is 1^2=1. Area of the n-th fractal is the area of the previous one + the areas of the newly arisen squares. The number of these new squares is equal to N(n-1), since there is one new square at every side of the previous fractal. The area of one such new square is L(n)^2, since the length of one side of this square is L(n).

So this

A(1)=1:
A(n)=A(n-1)+N(n-1)*L(n)^2:

can be the answer. But this is a recurrence which can be solved, so we can get A(n) as a function of n which is more practical. That's why I used the ?rsolve command just to get a different (more practical) formula for the same thing.

Thank you. That's nice and also faster than the previous way.

Thank you. That's nice and also faster than the previous way.

@Axel Vogt 

Thanks. The problem of those zeros is that it was not possible to measure some values there because they were probably too low, but it's important.

@Axel Vogt 

Thanks. The problem of those zeros is that it was not possible to measure some values there because they were probably too low, but it's important.

Thanks.

Here is the link (for others) to DirectSearch package v.2: http://www.maplesoft.com/applications/view.aspx?SID=101333

Thanks.

Here is the link (for others) to DirectSearch package v.2: http://www.maplesoft.com/applications/view.aspx?SID=101333

Do you need to use those `#msub(mi("n"),mo("&uminus0;"))`? I changed it to n0, n1, n2 and got one step further. :)

pde1 := diff(n1(x, t), t) = diff(n1(x, t), x)+n0(x, t)-n1(x, t);
pde2 := diff(n0(x, t), t) = 0.3139e-2*(diff(n0(x, t), x, x))-2*n0(x, t)+n2(x, t)+n1(x, t);
pde3 := diff(n2(x, t), t) = -(diff(n2(x, t), x))+n0(x, t)-n2(x, t);

Then I got the error: Error, (in pdsolve/numeric) initial/boundary conditions must be defined at one or two points for each independent variable

I think that the solver needs the boundary conditions for each variable at the same point(s). When I change the ibc to

ibc := n1(x, 0) = 0, n0(x, 0) = 0, n2(x, 0) = 0, n1(10000000, t) = 0, n0(1000000, t) = 0, n2(1000000, t) = 1, (D[1](n0))(1000000, t) = 0;

then

sol := pdsolve({pde1, pde2, pde3}, {ibc}, numeric);

works.