I encountered a memory problem in Mathematica, when I tried to process my experimental data. It's a system of two differential equations and I need to find most suitable parameters. Unfortunately I am not a Pro in Mathematica, so the program used a lot of memory, when the parameter epsilon is more than 0.4. When it less than 0.4, the program work properly. The command 'historylength = 0' and attempts to reduce the Accuracy Goal and WorkingPrecision didn`t help. I can't use ' clear Cache ', because there isnt a circle. I'm trying to understand what mistakes I made, and how I may limit the memory usage. I have already bought extra-RAM, now its 4GB, and now I haven't free memory-slots in motherboard

 Remove["Global`*"]; 
    T=13200; 
    L = 0.085; 
    e = 0.41; 
    v = 0.000557197; 
    q = 0.1618; 
    C0 = 0.0256; 
    R = 0.00075; 

    data = {{L,600,0.141124587},{L,1200,0.254134509},{L,1800,0.342888644},
            {L,2400,0.424476295},{L,3600,0.562844542},{L,4800,0.657111356},
            {L,6000,0.75137817},{L,7200,0.815876516},{L,8430,0.879823594},
            {L,9000,0.900771775},{L,13200,1}};

    model[(De_)?NumberQ, (Kf_)?NumberQ, (Y_)?NumberQ] := 
       model[De, Kf, Y] =  yeld /.Last[Last[
         NDSolve[{
          v (Ci^(1,0))[z,t]+(Ci^(0,1))[z,t]==
                      -((3 (1-e) Kf (Ci[z,t]-C0))/
                         (R e (1-(R Kf (1-R/r[z,t]))/De))),
         (r^(0,1))[z,t]==
                      (R^2 Kf (Ci[z,t]-C0))/
                         (q r[z,t]^2 (1-(R Kf (1-R/r[z,t]))/De)),
         (yeld^(0,1))[z,t]==
                      Y*(v e Ci[z,t])/(L q (1-e)),
    r[z,0]==R,
    Ci[z,0]==0,
    Ci[0,t]==0,
    yeld[z,0]==0},
    {r[z,t],Ci[z,t],yeld},{z,0,L},{t,0,T}]]]


    fit = FindFit[data, 
          {model[De, Kf, Y][z, t],
          {Y > 0.97, Y < 1.03, Kf > 10^-6, Kf < 10^-4, De > 10^-13, De < 10^-9}},       
          {{De,7*10^-13}, {Kf,  10^-5}, {Y, 1}},
          {z, t}, Method -> NMinimize]

   data = {{600,0.141124587},{1200,0.254134509},{1800,0.342888644},
        {2400,0.424476295},{3600,0.562844542},{4800,0.657111356},
        {6000,0.75137817},{7200,0.815876516},{8430,0.879823594},
        {9000,0.900771775},{13200,1}};

    YYY = model[
                De /. fit[[1]], 
                Kf /. fit[[2]],
                 Y /. fit[[3]]]; 

    Show[Plot[Evaluate[YYY[L,t]],{t,0,T},PlotRange->All],
    ListPlot[data,PlotStyle->Directive[PointSize[Medium],Red]]]

the link on the .nb file

http://www.4shared.com/folder/249TSjlz/_online.html