differential_equations plugged.
[scilab.git] / scilab / modules / differential_equations / sci_gateway / cpp / sci_ode.cpp
1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2011 - DIGITEO - Cedric DELAMARRE
4 *
5 * This file must be used under the terms of the CeCILL.
6 * This source file is licensed as described in the file COPYING, which
7 * you should have received as part of this distribution.  The terms
8 * are also available at
9 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10 *
11 */
12 /*--------------------------------------------------------------------------*/
13
14 #include "differential_equations_gw.hxx"
15 #include "function.hxx"
16 #include "double.hxx"
17 #include "string.hxx"
18 #include "list.hxx"
19 #include "callable.hxx"
20 #include "differentialequationfunctions.hxx"
21 #include "runvisitor.hxx"
22
23 extern "C"
24 {
25 #include "localization.h"
26 #include "Scierror.h"
27 #include "scifunctions.h"
28 #include "elem_common.h"
29 #include "sci_warning.h"
30 #include "sciprint.h"
31 #include "common_structure.h"
32 #include "checkodeerror.h"
33 }
34 /*--------------------------------------------------------------------------*/
35
36 types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, types::typed_list &out)
37 {
38     // Methode
39     types::String* pStrType     = NULL;
40     wchar_t* wcsType            = L"lsoda";
41     int meth                    = 0;
42
43     // y0
44 //    types::Polynom* pPolyY0     = NULL;
45     types::Double* pDblY0       = NULL;
46     double* pdYData             = NULL; // contain y0 following by all args data in list case.
47     int sizeOfpdYData           = 0;
48
49     // Other input args
50     types::Double* pDblT0       = NULL;
51     types::Double* pDblT        = NULL;
52     types::Double* pDblRtol     = NULL;
53     types::Double* pDblAtol     = NULL;
54     types::Double* pDblNg       = NULL;
55     types::Double* pDblW        = NULL;
56     types::Double* pDblIw       = NULL;
57
58     // %ODEOPTIONS
59     types::Double* pDblOdeOptions = NULL;
60
61     // Result
62     types::Double* pDblYOut     = NULL;
63 //    types::Double* pDblTOut     = NULL;
64 //    types::Polynom* pPolyYOut   = NULL;
65
66     // Indicate if the function is given.
67     bool bFuncF     = false;
68     bool bFuncJac   = false;
69     bool bFuncG     = false;
70
71     int iPos        = 0; // Position in types::typed_list in
72
73     int sizeOfYSize = 1;
74     int* YSize      = NULL;    // YSize(1) = size of y0, 
75                                // YSize(n) = size of Args(n) in list case.
76
77     C2F(eh0001).mesflg = 1; // flag to control printing of error messages in lapack routine.
78                             // 1 means print, 0 means no printing.
79     C2F(eh0001).lunit = 6;  // 6 = stdout
80
81     int one = 1; // used in C2F(dcopy)
82
83     // For root methode
84     int* jroot = NULL;
85
86 // *** check the minimal number of input args. ***
87     if(in.size() < 4)
88     {
89         ScierrorW(77, _W("%ls: Wrong number of input argument(s): %d expected.\n"), L"ode", 4);
90         return types::Function::Error;
91     }
92
93 // *** Get the methode. ***
94     if(in[0]->isString())
95     {
96         pStrType = in[0]->getAs<types::String>();
97         wcsType = pStrType->get(0);
98         iPos++;
99     }
100
101     if(iPos)
102     {
103         if(wcscmp(wcsType, L"adams") == 0)
104         {
105             meth = 1;
106         }
107         else if(wcscmp(wcsType, L"stiff") == 0)
108         {
109             meth = 2;
110         }
111         else if(wcscmp(wcsType, L"root") == 0)
112         {
113             meth = 3;
114         }
115         else if(wcscmp(wcsType, L"discrete") == 0)
116         {
117             meth = 4;
118         }
119         else if(wcscmp(wcsType, L"rk") == 0)
120         {
121             meth = 5;
122         }
123         else if(wcscmp(wcsType, L"rkf") == 0)
124         {
125             meth = 6;
126         }
127         else if(wcscmp(wcsType, L"fix") == 0)
128         {
129             meth = 7;
130         }
131         else
132         {
133             ScierrorW(999, _W("%ls: Wrong value for input argument #%d : It must be one of the following strings : adams, stiff, rk, rkf, fix, root or discrete.\n"), L"ode", 1);
134             return types::Function::Error;
135         }
136     }
137
138 // *** check number of output args according the methode. ***
139     if(meth < 3)
140     {
141         if(_iRetCount != 1 && _iRetCount != 3)
142         {
143             ScierrorW(78, _W("%ls: Wrong number of output argument(s): %d or %d expected.\n"), L"ode", 1, 3);
144             return types::Function::Error;
145         }
146     }
147     else if(meth == 3)
148     {
149         if(_iRetCount == 3 || _iRetCount > 4)
150         {
151             ScierrorW(78, _W("%ls: Wrong number of output argument(s): %d, %d or %d expected.\n"), L"ode", 1, 2, 4);
152             return types::Function::Error;
153         }
154     }
155     else // meth > 3
156     {
157         if(_iRetCount != 1)
158         {
159             ScierrorW(78, _W("%ls: Wrong number of output argument(s): %d expected.\n"), L"ode", 1);
160             return types::Function::Error;
161         }
162     }
163
164 // *** check type of input args and get it. ***
165     // y0
166     if(in[iPos]->isDouble())
167     {
168         pDblY0 = in[iPos]->getAs<types::Double>();
169         if(pDblY0->isComplex())
170         {
171             ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A real matrix expected.\n"), L"ode", iPos+1);
172             return types::Function::Error;
173         }
174
175         if(pDblY0->getCols() != 1)
176         {
177             ScierrorW(999, _W("%ls: Wrong size for input argument #%d : A real colunm vector expected (n x 1).\n"), L"ode", iPos+1);
178             return types::Function::Error;
179         }
180     }
181 /*
182     else if(in[iPos]->isPoly())
183     {
184         pPolyY0 = in[iPos]->getAs<types::Polynom>();
185         if(pPolyY0->isComplex())
186         {
187             ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A real polynom expected.\n"), L"ode", iPos+1);
188             return types::Function::Error;
189         }
190
191         if(pPolyY0->isScalar() == false)
192         {
193             ScierrorW(999, _W("%ls: Wrong size for input argument #%d : A real scalar polynom expected.\n"), L"ode", iPos+1);
194             return types::Function::Error;
195         }
196
197         double* dbl = (double*)malloc(pPolyY0->getCoef()->getSize() * sizeof(double));
198         vTransposeRealMatrix(   pPolyY0->getCoef()->get(),
199                                 pPolyY0->getCoef()->getRows(),
200                                 pPolyY0->getCoef()->getCols(),
201                                 dbl);
202
203         pDblY0 = new types::Double(pPolyY0->getCoef()->getCols(), pPolyY0->getCoef()->getRows());
204         pDblY0->set(dbl);
205     }
206 */
207     else
208     {
209         ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A matrix expected.\n"), L"ode", iPos+1);
210         return types::Function::Error;
211     }
212
213     // t0
214     iPos++;
215     if(in[iPos]->isDouble() == false)
216     {
217         ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A scalar expected.\n"), L"ode", iPos+1);
218         return types::Function::Error;
219     }
220
221     pDblT0 = in[iPos]->getAs<types::Double>();
222
223     if(pDblT0->isScalar() == false)
224     {
225         ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A scalar expected.\n"), L"ode", iPos+1);
226         return types::Function::Error;
227     }
228
229     // t
230     iPos++;
231     if(in[iPos]->isDouble() == false)
232     {
233         ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A matrix expected.\n"), L"ode", iPos+1);
234         return types::Function::Error;
235     }
236
237     pDblT = in[iPos]->getAs<types::Double>();
238
239     // get next inputs
240     DifferentialEquationFunctions* deFunctionsManager = new DifferentialEquationFunctions(L"ode");
241     DifferentialEquation::addDifferentialEquationFunctions(deFunctionsManager);
242     deFunctionsManager->setOdeYRows(pDblY0->getRows());
243     deFunctionsManager->setOdeYCols(pDblY0->getCols());
244
245     YSize = (int*)malloc(sizeOfYSize * sizeof(int));
246     *YSize = pDblY0->getSize();
247     pdYData = (double*)malloc(pDblY0->getSize() * sizeof(double));
248     C2F(dcopy)(YSize, pDblY0->get(), &one, pdYData, &one);
249
250     if(meth == 4)
251     {
252         if(in.size() != 5)
253         {
254             ScierrorW(77, _W("%ls: Wrong number of input argument(s): %d expected.\n"), L"ode", 5);
255             DifferentialEquation::removeDifferentialEquationFunctions();
256             free(pdYData);free(YSize);
257             return types::Function::Error;
258         }
259
260         if(in[4]->isCallable() == false && in[4]->isString() == false && in[4]->isList() == false)
261         {
262                 ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A function expected.\n"), L"ode", 5);
263                 DifferentialEquation::removeDifferentialEquationFunctions();
264                 free(pdYData);free(YSize);
265                 return types::Function::Error;
266         }
267     }
268
269     for(iPos++; iPos < in.size(); iPos++)
270     {
271         if(in[iPos]->isDouble())
272         {
273             if(pDblRtol == NULL && bFuncF == false)
274             {
275                 pDblRtol = in[iPos]->getAs<types::Double>();
276                 if(pDblRtol->getSize() != pDblY0->getSize() && pDblRtol->isScalar() == false)
277                 {
278                     ScierrorW(267,_W("%ls: Arg %d and arg %d must have equal dimensions.\n"), L"ode", pStrType ? 2 : 1, iPos+1);
279                     DifferentialEquation::removeDifferentialEquationFunctions();
280                     free(pdYData);free(YSize);
281                     return types::Function::Error;
282                 }
283             }
284             else if(pDblAtol == NULL && bFuncF == false)
285             {
286                 pDblAtol = in[iPos]->getAs<types::Double>();
287                 if(pDblAtol->getSize() != pDblY0->getSize() && pDblAtol->isScalar() == false)
288                 {
289                     ScierrorW(267,_W("%ls: Arg %d and arg %d must have equal dimensions.\n"), L"ode", pStrType ? 2 : 1, iPos+1);
290                     DifferentialEquation::removeDifferentialEquationFunctions();
291                     free(pdYData);free(YSize);
292                     return types::Function::Error;
293                 }
294             }
295             else if(pDblNg == NULL && bFuncF == true && meth == 3)
296             {
297                 pDblNg = in[iPos]->getAs<types::Double>();
298             }
299             else if(pDblW == NULL && bFuncF == true && (bFuncG == true || meth != 3))
300             {
301                 if(in.size() == iPos + 2)
302                 {
303                     if(in[iPos+1]->isDouble() == false)
304                     {
305                         ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A matrix expected.\n"), L"ode", iPos+2);
306                         DifferentialEquation::removeDifferentialEquationFunctions();
307                         free(pdYData);free(YSize);
308                         return types::Function::Error;
309                     }
310
311                     pDblW = in[iPos]->getAs<types::Double>();
312                     pDblIw = in[iPos+1]->getAs<types::Double>();
313                     iPos++;
314                 }
315                 else
316                 {
317                     ScierrorW(77, _W("%ls: Wrong number of input argument(s): %d expected.\n"), L"ode", iPos+2);
318                     DifferentialEquation::removeDifferentialEquationFunctions();
319                     free(pdYData);free(YSize);
320                     return types::Function::Error;
321                 }
322             }
323             else
324             {
325                 ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A function expected.\n"), L"ode", iPos+1);
326                 DifferentialEquation::removeDifferentialEquationFunctions();
327                 free(pdYData);free(YSize);
328                 return types::Function::Error;
329             }
330         }
331         else if(in[iPos]->isCallable())
332         {
333             types::Callable* pCall = in[iPos]->getAs<types::Callable>();
334             if(bFuncF == false)
335             {
336                 deFunctionsManager->setFFunction(pCall);
337                 bFuncF = true;
338             }
339             else if(bFuncJac == false && (pDblNg == NULL || meth != 3))
340             {
341                 deFunctionsManager->setJacFunction(pCall);
342                 bFuncJac = true;
343             }
344             else if(bFuncG == false && meth == 3)
345             {
346                 deFunctionsManager->setGFunction(pCall);
347                 bFuncG = true;
348             }
349             else
350             {
351                 ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A matrix expected.\n"), L"ode", iPos+1);
352                 DifferentialEquation::removeDifferentialEquationFunctions();
353                 free(pdYData);free(YSize);
354                 return types::Function::Error;
355             }
356         }
357         else if(in[iPos]->isString())
358         {
359             types::String* pStr = in[iPos]->getAs<types::String>();
360             bool bOK = false;
361
362             if(bFuncF == false)
363             {
364                 bOK = deFunctionsManager->setFFunction(pStr);
365                 bFuncF = true;
366             }
367             else if(bFuncJac == false && (pDblNg == NULL || meth != 3))
368             {
369                 bOK = deFunctionsManager->setJacFunction(pStr);
370                 bFuncJac = true;
371             }
372             else if(bFuncG == false && meth == 3)
373             {
374                 bOK = deFunctionsManager->setGFunction(pStr);
375                 bFuncG = true;
376             }
377             else
378             {
379                 ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A matrix expected.\n"), L"ode", iPos+1);
380                 DifferentialEquation::removeDifferentialEquationFunctions();
381                 free(pdYData);free(YSize);
382                 return types::Function::Error;
383             }
384
385             if(bOK == false)
386             {
387                 ScierrorW(50,_W("%ls: Subroutine not found: %ls\n"), L"ode", pStr->get(0));
388                 DifferentialEquation::removeDifferentialEquationFunctions();
389                 free(pdYData);free(YSize);
390                 return types::Function::Error;
391             }
392         }
393         else if(in[iPos]->isList())
394         {
395             types::List* pList = in[iPos]->getAs<types::List>();
396
397             if(pList->getSize() == 0)
398             {
399                 ScierrorW(50,_W("%ls: Argument #%d : Subroutine not found in list: %ls\n"), L"ode", iPos+1, L"(string empty)");
400                 DifferentialEquation::removeDifferentialEquationFunctions();
401                 free(pdYData);free(YSize);
402                 return types::Function::Error;
403             }
404
405             if(bFuncF && (bFuncJac || pDblNg) && (bFuncG || meth != 3))
406             {
407                 ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A matrix expected.\n"), L"ode", iPos+1);
408                 DifferentialEquation::removeDifferentialEquationFunctions();
409                 free(pdYData);free(YSize);
410                 return types::Function::Error;
411             }
412
413             if(pList->get(0)->isString())
414             {
415                 types::String* pStr = pList->get(0)->getAs<types::String>();
416                 bool bOK = false;
417
418                 if(bFuncF == false)
419                 {
420                     bFuncF = true;
421                     bOK = deFunctionsManager->setFFunction(pStr);
422                     sizeOfpdYData = *YSize;
423                 }
424                 else if(bFuncJac == false && (pDblNg == NULL || meth != 3))
425                 {
426                     bFuncJac = true;
427                     bOK = deFunctionsManager->setJacFunction(pStr);
428                     if(sizeOfpdYData == 0)
429                     {
430                         sizeOfpdYData = *YSize;
431                     }
432                 }
433                 else if(bFuncG == false && meth == 3)
434                 {
435                     bFuncG = true;
436                     bOK = deFunctionsManager->setGFunction(pStr);
437                     if(sizeOfpdYData == 0)
438                     {
439                         sizeOfpdYData = *YSize;
440                     }
441                 }
442
443                 if(bOK == false)
444                 {
445                     ScierrorW(50,_W("%ls: Argument #%d : Subroutine not found in list: %ls\n"), L"ode", iPos+1, pStr->get(0));
446                     DifferentialEquation::removeDifferentialEquationFunctions();
447                     free(pdYData);free(YSize);
448                     return types::Function::Error;
449                 }
450
451                 int* sizeTemp = YSize;
452                 int totalSize = sizeOfpdYData;
453
454                 YSize = (int*)malloc((sizeOfYSize + pList->getSize() - 1) * sizeof(int));
455                 memcpy(YSize, sizeTemp, sizeOfYSize * sizeof(int));
456
457                 std::vector<types::Double*> vpDbl;
458                 for(int iter = 0; iter < pList->getSize() - 1; iter++)
459                 {
460                     if(pList->get(iter + 1)->isDouble() == false)
461                     {
462                         ScierrorW(999, _W("%ls: Wrong type for input argument #%d : Argument %d in the list must be a matrix.\n"), L"ode", iPos+1, iter+1);
463                         DifferentialEquation::removeDifferentialEquationFunctions();
464                         free(pdYData);free(YSize);
465                         return types::Function::Error;
466                     }
467
468                     vpDbl.push_back(pList->get(iter+1)->getAs<types::Double>());
469                     YSize[sizeOfYSize + iter] = vpDbl[iter]->getSize();
470                     totalSize += YSize[sizeOfYSize + iter];
471                 }
472
473                 double* pdYDataTemp = pdYData;
474                 pdYData = (double*)malloc(totalSize * sizeof(double));
475                 C2F(dcopy)(&sizeOfpdYData, pdYDataTemp, &one, pdYData, &one);
476
477                 int position = sizeOfpdYData;
478                 for(int iter = 0; iter < pList->getSize()-1; iter++)
479                 {
480                     C2F(dcopy)(&YSize[sizeOfYSize + iter], vpDbl[iter]->get(), &one, &pdYData[position], &one);
481                     position += vpDbl[iter]->getSize();
482                 }
483                 vpDbl.clear();
484                 sizeOfpdYData = totalSize;
485                 sizeOfYSize += pList->getSize() - 1;
486                 free(pdYDataTemp);
487                 free(sizeTemp);
488             }
489             else if(pList->get(0)->isCallable())
490             {
491                 if(bFuncF == false)
492                 {
493                     bFuncF = true;
494                     deFunctionsManager->setFFunction(pList->get(0)->getAs<types::Callable>());
495                     for(int iter = 1; iter < pList->getSize(); iter++)
496                     {
497                         deFunctionsManager->setFArgs(pList->get(iter)->getAs<types::InternalType>());
498                     }
499                 }
500                 else if(bFuncJac == false && (pDblNg == NULL || meth != 3))
501                 {
502                     bFuncJac = true;
503                     deFunctionsManager->setJacFunction(pList->get(0)->getAs<types::Callable>());
504                     for(int iter = 1; iter < pList->getSize(); iter++)
505                     {
506                         deFunctionsManager->setJacArgs(pList->get(iter)->getAs<types::InternalType>());
507                     }
508                 }
509                 else if(bFuncG == false && meth == 3)
510                 {
511                     bFuncG = true;
512                     deFunctionsManager->setGFunction(pList->get(0)->getAs<types::Callable>());
513                     for(int iter = 1; iter < pList->getSize(); iter++)
514                     {
515                         deFunctionsManager->setGArgs(pList->get(iter)->getAs<types::InternalType>());
516                     }
517                 }
518             }
519             else
520             {
521                 ScierrorW(999, _W("%ls: Wrong type for input argument #%d : The first argument in the list must be a string or a function.\n"), L"ode", iPos+1);
522                 DifferentialEquation::removeDifferentialEquationFunctions();
523                 free(pdYData);free(YSize);
524                 return types::Function::Error;
525             }
526         }
527         else
528         {
529             ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A matrix or a function expected.\n"), L"ode", iPos+1);
530             DifferentialEquation::removeDifferentialEquationFunctions();
531             free(pdYData);free(YSize);
532             return types::Function::Error;
533         }
534     }
535
536     if(bFuncF == false)
537     {
538         int val = (meth == 3) ? 3 : 1;
539         ScierrorW(77, _W("%ls: Wrong number of input argument(s): %d expected.\n"), L"ode", in.size() + val);
540         DifferentialEquation::removeDifferentialEquationFunctions();
541         free(pdYData);free(YSize);
542         return types::Function::Error;
543     }
544     if(pDblNg == NULL && meth == 3)
545     {
546         ScierrorW(77, _W("%ls: Wrong number of input argument(s): %d expected.\n"), L"ode", in.size() + 2);
547         DifferentialEquation::removeDifferentialEquationFunctions();
548         free(pdYData);free(YSize);
549         return types::Function::Error;
550     }
551     if(bFuncG == false && meth == 3)
552     {
553         ScierrorW(77, _W("%ls: Wrong number of input argument(s): %d expected.\n"), L"ode", in.size() + 1);
554         DifferentialEquation::removeDifferentialEquationFunctions();
555         free(pdYData);free(YSize);
556         return types::Function::Error;
557     }
558
559 // *** Initialization. ***
560
561     int itol        = 1;
562     int itask       = 1; // itask = 1 for normal computation of output values of y at t = tout.
563     int istate      = 1; // istate = integer flag (input and output).  set istate = 1.
564     int iopt        = 0; // iopt   = 0 to indicate no optional inputs used.
565     int jt          = 2; // jacobian type indicator.  set jt = 2.
566     int ml          = -1;
567     int mu          = -1;
568
569     // work tab
570     double* rwork   = NULL;
571     int* iwork      = NULL;
572     int rworkSize   = 0;
573     int iworkSize   = 0;
574
575     // contain ls0001, lsa001 and eh0001 structures
576     double* dStructTab  = NULL;
577     int* iStructTab     = NULL;
578     int dStructTabSize  = 0;
579     int iStructTabSize  = 0;
580
581     int rwSize  = 0; // rwSize = dStructTab + rworkSize
582     int iwSize  = 0; // iwSize = iStructTab + iworkSize
583
584     // structures used by lsoda and lsode
585     double* ls0001d = &(C2F(ls0001).tret);
586     int* ls0001i    = &(C2F(ls0001).illin);
587     double* lsa001d = &(C2F(lsa001).tsw);
588     int* lsa001i    = &(C2F(lsa001).insufr);
589     double* lsr001d = &(C2F(lsr001).rownr3[0]);
590     int* lsr001i    = &(C2F(lsr001).lg0);
591     int* eh0001i    = &(C2F(eh0001).mesflg);
592
593     // get %ODEOPTIONS
594     types::InternalType* pIT = symbol::Context::getInstance()->get(symbol::Symbol(L"%ODEOPTIONS"));
595     if(pIT != NULL && pIT->isDouble())
596     {
597         pDblOdeOptions = pIT->getAs<types::Double>();
598         if(pDblOdeOptions->getSize() == 12)
599         {
600             iopt    = 1;
601             itask   = (int)pDblOdeOptions->get(0);
602             jt      = (int)pDblOdeOptions->get(5);
603             ml      = (int)pDblOdeOptions->get(10);
604             mu      = (int)pDblOdeOptions->get(11);
605         }
606     }
607
608     if(iopt == 1 && (pDblOdeOptions->get(4) > pDblOdeOptions->get(3))) // hmin > hmax ?
609     {
610         ScierrorW(9999, _W("%ls: Wrong value of hmin and hmax : hmin = %d is greater than hmax = %d.\n"), L"ode", pDblOdeOptions->get(4), pDblOdeOptions->get(3));
611         DifferentialEquation::removeDifferentialEquationFunctions();
612         free(pdYData);free(YSize);
613         return types::Function::Error;
614     }
615
616     if(jt < 0 || jt > 5)
617     {
618         ScierrorW(9999, _W("%ls: Wrong value of Jacobian type : A number between %d and %d expected.\n"), L"ode", 0, 5);
619         DifferentialEquation::removeDifferentialEquationFunctions();
620         free(pdYData);free(YSize);
621         return types::Function::Error;
622     }
623
624     if(iopt == 0 && bFuncJac)
625     {
626         jt = 1;
627     }
628
629     if(bFuncJac && (jt == 2 || jt == 5) && getWarningMode())
630     {
631         sciprintW(_W("%ls: Warning: Jacobian is given, but not used.\n"), L"ode");
632     }
633
634     if(bFuncJac == false && (jt == 1 || jt == 4))
635     {
636         if(getWarningMode())
637         {
638             sciprintW(_W("%ls: Warning: No Jacobian external given, but one is required by %ODEOPTIONS(6) value. Jacobian will be estimated.\n"), L"ode");
639         }
640
641         jt = 2;
642     }
643
644     //compute itol and set the tolerances rtol and atol.
645     double* rtol = NULL;
646     double* atol = NULL;
647
648     if(pDblRtol)
649     {
650         if(pDblRtol->isScalar())
651         {
652             rtol = (double*)malloc(sizeof(double));
653             *rtol = pDblRtol->get(0);
654         }
655         else
656         {
657             rtol = pDblRtol->get();
658             itol += 2;
659         }
660     }
661     else
662     {
663         rtol = (double*)malloc(sizeof(double));
664         if(meth == 6 || meth == 7)
665         {
666             *rtol = 1.e-3;
667         }
668         else
669         {
670             *rtol = 1.0e-7;
671         }
672     }
673
674     if(pDblAtol)
675     {
676         if(pDblAtol->isScalar())
677         {
678             atol = (double*)malloc(sizeof(double));
679             *atol = pDblAtol->get(0);
680         }
681         else
682         {
683             atol = pDblAtol->get();
684             itol ++;
685         }
686     }
687     else
688     {
689         atol = (double*)malloc(sizeof(double));
690         if(meth == 6 || meth == 7)
691         {
692             *atol = 1.e-4;
693         }
694         else
695         {
696             *atol = 1.0e-9;
697         }
698     }
699
700     // Compute rwork, iwork size.
701     // Create them.
702     int nyh = (*YSize);
703     if(pDblW) // structure ls0001 have been restored.
704     {
705         nyh = C2F(ls0001).nyh;
706     }
707
708     int mxordn = 12;
709     int mxords = 5;
710
711     if(iopt == 1)
712     {
713         mxordn = (int)pDblOdeOptions->get(7);
714         mxords = (int)pDblOdeOptions->get(8);
715     }
716
717     if(mxordn > 12 || mxords > 5 || mxordn < 1 || mxords < 1)
718     {
719         if(getWarningMode())
720         {
721             sciprintW(_W("%ls: Warning: Wrong value for maximun stiff/non-stiff order allowed :\nAt most %d for mxordn, %d for mxords and no null value for both expected.\nWrong value will be reduced to the default value.\n"), L"ode", 12, 5);
722         }
723
724         mxordn = 12;
725         mxords = 5;
726     }
727
728     int maxord  = mxords;
729     int lmat    = 0;
730     int lrn     = 0;
731     int lrs     = 0;
732
733     switch(meth)
734     {
735         case 3 : // lsodar (root)
736         {
737 //             lrn = 20 + nyh*(mxordn+1) + 3*neq + 3*ng,
738 //             lrs = 20 + nyh*(mxords+1) + 3*neq + lmat + 3*ng,
739 //          where
740 //             nyh    = the initial value of neq,
741 //             mxordn = 12, unless a smaller value is given as an
742 //                      optional input,
743 //             mxords = 5, unless a smaller value is given as an
744 //                      optional input,
745 //             lmat   = length of matrix work space..
746 //             lmat   = neq**2 + 2              if jt = 1 or 2,
747 //             lmat   = (2*ml + mu + 1)*neq + 2 if jt = 4 or 5.
748
749             lrn = 3 * (int)pDblNg->get(0);
750             lrs = lrn;
751             dStructTabSize = 246 - 241;
752             iStructTabSize = 59 - 50;
753             jroot = (int*)malloc((int)pDblNg->get(0) * sizeof(int));
754         }
755         case 0 : // lsoda
756         {
757 //             lrn = 20 + nyh*(mxordn+1) + 3*neq,
758 //             lrs = 20 + nyh*(mxords+1) + 3*neq + lmat,
759 //          where
760 //             nyh    = the initial value of neq,
761 //             mxordn = 12, unless a smaller value is given as an
762 //                      optional input,
763 //             mxords = 5, unless a smaller value is given as an
764 //                      optional input,
765 //             lmat   = length of matrix work space..
766 //             lmat   = neq**2 + 2              if jt = 1 or 2,
767 //             lmat   = (2*ml + mu + 1)*neq + 2 if jt = 4 or 5.
768
769             if(jt == 1 || jt == 2)
770             {
771                 lmat = (*YSize) * (*YSize) + 2;       //  if jt = 1 or 2,
772             }
773             else if(jt == 4 || jt == 5)
774             {   //ml and mu = -1 in all cases
775                 lmat = (2 * ml + mu + 1) * (*YSize) + 2; //  if jt = 4 or 5. 
776             }
777
778             lrn += 20 + nyh * (mxordn + 1) + 3 * (*YSize);
779             lrs += 20 + nyh * (mxords + 1) + 3 * (*YSize) + lmat;
780
781             rworkSize   = max(lrn, lrs);
782             iworkSize   = 20 + *YSize;     
783
784             dStructTabSize += 241;
785             iStructTabSize += 50;
786
787             break;
788         }
789         case 1 : // lsode (adams) (non stiff)
790         {
791             maxord = mxordn;
792         }
793         case 2 : // lsode (stiff)
794         {
795 //          20 + nyh*(maxord + 1) + 3*neq + lmat
796 //        where
797 //          nyh    = the initial value of neq,
798 //          maxord = 12 (if meth = 1) or 5 (if meth = 2) (unless a
799 //                   smaller value is given as an optional input),
800 //          lmat   = 0             if miter = 0,
801 //          lmat   = neq**2 + 2    if miter is 1 or 2,
802 //          lmat   = neq + 2       if miter = 3, and
803 //          lmat   = (2*ml+mu+1)*neq + 2 if miter is 4 or 5.
804
805             if(jt == 1 || jt == 2)
806             {
807                 lmat = (*YSize) * (*YSize) + 2;
808             }
809             else if(jt == 3)
810             {
811                 lmat = (*YSize) + 2;
812             }
813             else if(jt == 4 || jt == 5)
814             {
815                 lmat = (2 * ml + mu + 1) * (*YSize) + 2;
816             }
817
818             rworkSize = 20 + nyh * (maxord + 1) + 3 * (*YSize) + lmat;
819             iworkSize = 20; // if jt = 0 or 3
820             if(jt == 1 || jt == 2 || jt == 4 || jt == 5)// iSize = 20 + neq  otherwise
821             {
822                 iworkSize += (*YSize);
823             }
824
825             dStructTabSize = 219;
826             iStructTabSize = 41;
827             break;
828         }
829         case 4 : // lsdisc (discrete)
830         {
831             rworkSize = *YSize;
832             iworkSize = 1;
833             break;
834         }
835         case 5 : // lsrgk (rk)
836         {
837             rworkSize = 3 * (*YSize);
838             iworkSize = 1;
839             break;
840         }
841         case 6 : // rkf45 (rkf)
842         {
843             rworkSize = 3 + 8 * (*YSize);
844             iworkSize = 5;
845             break;
846         }
847         case 7 : // rksimp (fix)
848         {
849             rworkSize = 3 + 8 * (*YSize);
850             iworkSize = 1;
851             break;
852         }
853     }
854
855     rwSize = rworkSize + dStructTabSize;
856     iwSize = iworkSize + iStructTabSize;
857
858     rwork = (double*)malloc(rworkSize * sizeof(double));
859     iwork = (int*)malloc(iworkSize * sizeof(int));
860
861     if(meth < 4)
862     {
863         if(pDblW && pDblIw)
864         {
865             if(pDblW->getSize() != rwSize || pDblIw->getSize() != iwSize)
866             {
867                 ScierrorW(9999, _W("%ls: Wrong size for w and iw : w = %d and iw = %d expected.\n"), L"ode", rwSize, iwSize);
868                 DifferentialEquation::removeDifferentialEquationFunctions();
869                 free(pdYData);free(YSize);
870                 free(rwork);free(iwork);
871                 if(jroot) free(jroot);
872                 if(itol == 1 || itol == 3) free(atol);
873                 if(itol < 3) free(rtol);
874                 return types::Function::Error;
875             }
876
877             istate = 2; // 1 means this is the first call | 2  means this is not the first call
878
879             dStructTab = (double*)malloc(dStructTabSize * sizeof(double));
880             iStructTab = (int*)malloc(iStructTabSize * sizeof(int));
881
882             C2F(dcopy)(&rworkSize, pDblW->get(), &one, rwork, &one);
883             C2F(dcopy)(&dStructTabSize, pDblW->get() + rworkSize, &one, dStructTab, &one);
884
885             for(int i = 0; i < iworkSize; i++)
886             {
887                 iwork[i] = (int)pDblIw->get(i);
888             }
889
890             for(int i = 0; i < iStructTabSize; i++)
891             {
892                 iStructTab[i] = (int)pDblIw->get(i + iworkSize);
893             }
894         }
895         else
896         {
897             //if iopt = 0 default value are used without set it.
898             if(iopt == 1)
899             {
900                 if(itask >= 4)
901                 {
902                     rwork[0] = pDblOdeOptions->get(1); // tcrit
903                 }
904
905                 rwork[4] = pDblOdeOptions->get(2);      // h0 =0
906                 rwork[5] = pDblOdeOptions->get(3);      // hmax = %inf
907                 rwork[6] = pDblOdeOptions->get(4);      // hmin = 0
908                 iwork[0] = (int)pDblOdeOptions->get(10);// ml = -1
909                 iwork[1] = (int)pDblOdeOptions->get(11);// mu = -1
910                 iwork[4] = (int)pDblOdeOptions->get(9); // ixpr = 0
911                 iwork[5] = (int)pDblOdeOptions->get(6); // mxstep = 500
912                 iwork[6] = 10;  // mxhnil = 10 maximum number of messages printed per problem
913                 iwork[7] = (int)pDblOdeOptions->get(7); // mxordn = 12
914                 iwork[8] = (int)pDblOdeOptions->get(8); // mxords = 5
915             }
916         }
917     }
918
919     if(pDblOdeOptions && pDblOdeOptions->get(9) == 1)
920     {
921         sciprintW(_W("itask = %d\tmeth = %d\tjactyp = %d\tml = %d\tmu = %d\tiopt = %d\n"), itask, meth, jt, ml, mu, iopt);
922         sciprintW(_W("tcrit = %lf\th0 = %lf\thmax = %lf\thmin = %lf\n"), pDblOdeOptions->get(1), pDblOdeOptions->get(2), pDblOdeOptions->get(3), pDblOdeOptions->get(4));
923     }
924
925     if(meth < 4)
926     {
927         if(pDblW && pDblIw)
928         {
929             int dPos    = 0;
930             int iPos    = 0;
931             int dSize   = 219;
932
933             //restore ls0001 and eh0001 from w (rwork) and iw (iwork).
934             C2F(dcopy)(&dSize, dStructTab, &one, ls0001d, &one);
935             memcpy(ls0001i, iStructTab, 39 * sizeof(int));
936
937             dPos = dSize;
938             iPos = 39;
939
940             //restore lsa001 from w (rwork) and iw (iwork).
941             if(meth == 0 || meth == 3)
942             {
943                 dSize = 22;
944                 C2F(dcopy)(&dSize, dStructTab + dPos, &one, lsa001d, &one);
945                 memcpy(lsa001i, &iStructTab[iPos], 9 * sizeof(int));
946
947                 dPos += dSize;
948                 iPos += 9;
949             }
950
951             //restore lsr001
952             if(meth == 3)
953             {
954                 dSize = 5;
955                 C2F(dcopy)(&dSize, dStructTab + dPos, &one, lsr001d, &one);
956                 memcpy(lsr001i, &iStructTab[iPos], 9 * sizeof(int));
957             }
958         }
959     }
960
961 // *** Perform operation. ***
962     double ret = 0;
963     double t0 = pDblT0->get(0);
964     bool bOneStep = false;
965
966     if(itask == 5 || itask == 3 || itask == 2)
967     {
968         bOneStep = true;
969         if(getWarningMode() && pDblT->isScalar() == false)
970         {
971             sciprintW(_W("itask = %d : At most one value of t is allowed, the last element of t is used.\n"), itask);
972         }
973     }
974
975     if(bOneStep)
976     {
977         std::list<double*> pDblYOutList = std::list<double*>();
978         std::list<double> pDblTOutList = std::list<double>();
979         int iLastT = pDblT->getSize() - 1;
980         double t = pDblT->get(iLastT);
981
982         do
983         {
984             wchar_t* strMeth;
985             switch(meth)
986             {
987                 case 0 : // lsoda
988                 {
989                     strMeth = _W("lsoda");
990                     ret = C2F(lsoda)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt);
991                     break;
992                 }
993                 case 1 : // lsode (adams)
994                 case 2 : // lsode (stiff)
995                 {
996                     strMeth = _W("lsode");
997                     int jacType = 10*meth + jt;
998                     ret = C2F(lsode)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jacType);
999                     break;
1000                 }
1001                 case 3 : // lsodar
1002                 {
1003                     strMeth = _W("lsodar");
1004                     int ng = (int)pDblNg->get(0);
1005                     ret = C2F(lsodar)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt, ode_g, &ng, jroot);
1006                     break;
1007                 }
1008                 case 4 : // lsdisc (discrete)
1009                 {
1010                     strMeth = _W("lsdisc");
1011                     ret = C2F(lsdisc)(ode_f, YSize, pdYData, &t0, &t, rwork, &rworkSize, &istate);
1012                     break;
1013                 }
1014                 case 5 : // lsrgk (rk)
1015                 {
1016                     strMeth = _W("lsrgk");
1017                     ret = C2F(lsrgk)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1018                     break;
1019                 }
1020                 case 6 : // rkf45 (rkf)
1021                 {
1022                     strMeth = _W("rkf45");
1023                     ret = C2F(rkf45)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1024                     break;
1025                 }
1026                 case 7 : // rksimp (fix)
1027                 {
1028                     strMeth = _W("rksimp");
1029                     ret = C2F(rksimp)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1030                     break;
1031                 }
1032             }
1033             // check error
1034             int err = checkOdeError(meth, istate);
1035             if(err == 1) // error case
1036             {
1037                 ScierrorW(999, _W("%ls: %ls exit with state %d.\n"), L"ode", strMeth, istate);
1038                 DifferentialEquation::removeDifferentialEquationFunctions();
1039                 free(pdYData);free(YSize);
1040                 free(rwork);free(iwork);
1041                 if(jroot) free(jroot);
1042                 if(dStructTab) free(dStructTab);
1043                 if(iStructTab) free(iStructTab);
1044                 if(itol == 1 || itol == 3) free(atol);
1045                 if(itol < 3) free(rtol);
1046                 return types::Function::Error;
1047             }
1048
1049             pDblYOutList.push_back(pdYData);
1050             pDblTOutList.push_back(t0);
1051
1052             if(err == 2) // warning case
1053             {
1054                 if(getWarningMode())
1055                 {
1056                     sciprintW(_W("Integration was stoped at t = %lf.\n"),t0);
1057                 }
1058                 break;
1059             }
1060
1061             if(meth == 3 && istate == 3 && getWarningMode())
1062             {
1063             // istate == 3  means the integration was successful, and one or more
1064             //              roots were found before satisfying the stop condition
1065             //              specified by itask.  see jroot.
1066
1067                 sciprintW(_W("%ls: Warning: At t = %lf, y is a root, jroot = "), L"ode", t0);
1068                 for(int k = 0; k < pDblNg->get(0); k++)
1069                 {
1070                     sciprintW(_W("\t%d"), jroot[k]);
1071                 }
1072                 sciprintW(_W("\n"));
1073                 break;
1074             }
1075         }
1076         while(t0 < t);
1077
1078         int iSizeList = (int)pDblYOutList.size();
1079 /*
1080         if(pPolyY0) // pPolyY0 is scalar
1081         {
1082             // result format = [2 x iSizeList]
1083             // first lime is the time of result stored in second line.
1084             int size = 2 * iSizeList;
1085             int* iRanks = (int*)malloc(size * sizeof(int));
1086
1087             for(int i = 0; i < iSizeList; i++)
1088             {
1089                 iRanks[i * 2] = 1; // time rank
1090                 iRanks[i * 2 + 1] = pPolyY0->getMaxRank(); // result rank
1091             }
1092
1093             pPolyYOut = new types::Polynom(pPolyY0->getVariableName(), 2, iSizeList, iRanks);
1094             pDblYOut = new types::Double(pDblY0->getRows(), 1);
1095             pDblTOut = new types::Double(1, 1);
1096
1097             for(int i = 0; i < iSizeList; i++)
1098             {
1099                 // pDblY0 contain pPolyY0 (ie : pPolyY0 = 2 + x => pDblY0 = [2 1])
1100                 for(int j = 0; j < pDblY0->getRows(); j++)
1101                 {
1102                     pDblYOut->set(j, pDblYOutList.front()[j]);
1103                 }
1104                 pDblTOut->set(0, pDblTOutList.front());
1105
1106                 pPolyYOut->setCoef(i * 2, pDblTOut);
1107                 pPolyYOut->setCoef(i * 2 + 1, pDblYOut);
1108
1109                 pDblYOutList.pop_front();
1110                 pDblTOutList.pop_front();
1111             }
1112         }
1113         else
1114         {
1115 */
1116             pDblYOut = new types::Double(pDblY0->getRows() + 1, iSizeList);
1117             for(int i = 0; i < iSizeList; i++)
1118             {
1119                 pDblYOut->set(i * (*YSize + 1), pDblTOutList.front());
1120                 for(int j = 0; j < *YSize; j++)
1121                 {
1122                     pDblYOut->set(i * (*YSize + 1) + (j + 1), pDblYOutList.front()[j]);
1123                 }
1124                 pDblYOutList.pop_front();
1125                 pDblTOutList.pop_front();
1126             }
1127 //        }
1128     }
1129     else
1130     {
1131         // Create result
1132 /*        if(pPolyY0)
1133         {
1134             int size = pDblT->getSize();
1135             int* iRanks = (int*)malloc(size * sizeof(int));
1136             for(int i = 0; i < size; i++)
1137             {
1138                 iRanks[i] = pPolyY0->getMaxRank();
1139             }
1140
1141             pPolyYOut = new types::Polynom(pPolyY0->getVariableName(), 1, pDblT->getSize(), iRanks);
1142
1143             pDblYOut = new types::Double(pDblY0->getRows(), 1);
1144         }
1145         else
1146         {
1147 */
1148             pDblYOut = new types::Double(pDblY0->getRows(), pDblT->getSize());
1149 //        }
1150         bool bBreak = false;
1151         for(int i = 0; i < pDblT->getSize(); i++)
1152         {
1153             double t = pDblT->get(i);
1154             wchar_t* strMeth;
1155
1156             if(itask >= 4 && t > rwork[0]) // rwork[0] => tcrit
1157             {
1158                 t = rwork[0];
1159                 bBreak = true;
1160             }
1161
1162             switch(meth)
1163             {
1164                 case 0 : // lsoda
1165                 {
1166                     strMeth = _W("lsoda");
1167                     ret = C2F(lsoda)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt);
1168                     break;
1169                 }
1170                 case 1 : // lsode (adams)
1171                 case 2 : // lsode (stiff)
1172                 {
1173                     strMeth = _W("lsode");
1174                     int jacType = 10*meth + jt;
1175                     ret = C2F(lsode)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jacType);
1176                     break;
1177                 }
1178                 case 3 : // lsodar
1179                 {
1180                     strMeth = _W("lsodar");
1181                     int ng = (int)pDblNg->get(0);
1182                     ret = C2F(lsodar)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt, ode_g, &ng, jroot);
1183                     break;
1184                 }
1185                 case 4 : // lsdisc (discrete)
1186                 {
1187                     strMeth = _W("lsdisc");
1188                     ret = C2F(lsdisc)(ode_f, YSize, pdYData, &t0, &t, rwork, &rworkSize, &istate);
1189                     break;
1190                 }
1191                 case 5 : // lsrgk (rk)
1192                 {
1193                     strMeth = _W("lsrgk");
1194                     ret = C2F(lsrgk)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1195                     break;
1196                 }
1197                 case 6 : // rkf45 (rkf)
1198                 {
1199                     strMeth = _W("rkf45");
1200                     ret = C2F(rkf45)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1201                     break;
1202                 }
1203                 case 7 : // rksimp (fix)
1204                 {
1205                     strMeth = _W("rksimp");
1206                     ret = C2F(rksimp)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1207                     break;
1208                 }
1209             }
1210             // check error
1211             int err = checkOdeError(meth, istate);
1212             if(err == 1) // error case
1213             {
1214                 ScierrorW(999, _W("%ls: %ls exit with state %d.\n"), L"ode", strMeth, istate);
1215                 DifferentialEquation::removeDifferentialEquationFunctions();
1216                 free(pdYData);free(YSize);
1217                 free(rwork);free(iwork);
1218                 if(jroot) free(jroot);
1219                 if(dStructTab) free(dStructTab);
1220                 if(iStructTab) free(iStructTab);
1221                 if(itol == 1 || itol == 3) free(atol);
1222                 if(itol < 3) free(rtol);
1223                 return types::Function::Error;
1224             }
1225
1226 /*            if(pPolyYOut)
1227             {
1228                 for(int j = 0; j < pDblY0->getRows(); j++)
1229                 {
1230                     pDblYOut->set(j, pdYData[j]);
1231                 }
1232                 pPolyYOut->setCoef(i, pDblYOut);
1233             }
1234             else
1235             {
1236 */
1237                 for(int j = 0; j < *YSize; j++)
1238                 {
1239                     pDblYOut->set(i*(*YSize) + j, pdYData[j]);
1240                 }
1241 //            }
1242
1243             if(err == 2) // warning case
1244             {
1245                 if(getWarningMode())
1246                 {
1247                     sciprintW(_W("Integration was stoped at t = %lf.\n"),t0);
1248                 }
1249
1250                 types::Double* pDblYOutTemp = pDblYOut;
1251                 pDblYOut = new types::Double(pDblYOutTemp->getRows(), i+1);
1252                 for(int k = 0; k < pDblYOut->getSize(); k++)
1253                 {
1254                     pDblYOut->set(k, pDblYOutTemp->get(k));
1255                 }
1256                 delete pDblYOutTemp;
1257
1258                 break;
1259             }
1260
1261             if(meth == 3 && istate == 3 && getWarningMode())
1262             {
1263             // istate == 3  means the integration was successful, and one or more
1264             //              roots were found before satisfying the stop condition
1265             //              specified by itask.  see jroot.
1266
1267                 sciprintW(_W("%ls: Warning: At t = %lf, y is a root, jroot = "), L"ode", t0);
1268                 for(int k = 0; k < pDblNg->get(0); k++)
1269                 {
1270                     sciprintW(_W("\t%d"), jroot[k]);
1271                 }
1272                 sciprintW(_W("\n"));
1273
1274                 types::Double* pDblYOutTemp = pDblYOut;
1275                 pDblYOut = new types::Double(pDblYOutTemp->getRows(), i+1);
1276                 for(int k = 0; k < pDblYOut->getSize(); k++)
1277                 {
1278                     pDblYOut->set(k, pDblYOutTemp->get(k));
1279                 }
1280                 delete pDblYOutTemp;
1281                 break;
1282             }
1283
1284             if(itask >= 4 && bBreak)
1285             {
1286                 types::Double* pDblYOutTemp = pDblYOut;
1287                 pDblYOut = new types::Double(pDblYOutTemp->getRows(), i+1);
1288                 for(int k = 0; k < pDblYOut->getSize(); k++)
1289                 {
1290                     pDblYOut->set(k, pDblYOutTemp->get(k));
1291                 }
1292                 delete pDblYOutTemp;
1293                 break;
1294             }
1295         }
1296     }
1297
1298     if(meth < 4)
1299     {
1300         if(_iRetCount > 2)//save ls0001 and eh0001 following w and iw.
1301         {
1302             int dPos    = 0;
1303             int iPos    = 0;
1304             int dSize   = 219;
1305
1306             if(dStructTab == NULL)
1307             {
1308                 dStructTab = (double*)malloc(dStructTabSize * sizeof(double));
1309                 iStructTab = (int*)malloc(iStructTabSize * sizeof(int));
1310             }
1311
1312             C2F(dcopy)(&dSize, ls0001d, &one, dStructTab, &one);
1313             memcpy(iStructTab, ls0001i, 39 * sizeof(int));
1314
1315             dPos = dSize;
1316             iPos = 39;
1317
1318             //save lsa001
1319             if(meth == 0 || meth == 3)
1320             {
1321                 dSize = 22;
1322                 C2F(dcopy)(&dSize, lsa001d, &one, dStructTab + dPos, &one);
1323                 memcpy(&iStructTab[iPos], lsa001i, 9 * sizeof(int));
1324
1325                 dPos += dSize;
1326                 iPos += 9;
1327             }
1328
1329             //save lsr001
1330             if(meth == 3)
1331             {
1332                 dSize = 5;
1333                 C2F(dcopy)(&dSize, lsr001d, &one, dStructTab + dPos, &one);
1334                 memcpy(&iStructTab[iPos], lsr001i, 9 * sizeof(int));
1335
1336                 iPos += 9;
1337             }
1338
1339             memcpy(&iStructTab[iPos], eh0001i, 2 * sizeof(int));
1340         }
1341     }
1342
1343 // *** Return result in Scilab. ***
1344 /*
1345     if(pPolyYOut)
1346     {
1347         out.push_back(pPolyYOut); // y
1348     }
1349     else
1350     {
1351 */
1352         out.push_back(pDblYOut); // y
1353 //    }
1354
1355     if(meth == 3 && _iRetCount >= 2)
1356     {
1357         int sizeOfRd = 1;
1358         int k = 0;
1359
1360         for(int i = 0; i < pDblNg->get(0); i++)
1361         {
1362             if(jroot[i])
1363             {
1364                 sizeOfRd++;
1365             }
1366         }
1367
1368         types::Double* pDblRd = new types::Double(1, sizeOfRd);
1369         //rd : The first entry contains the stopping time.
1370         pDblRd->set(0, C2F(lsr001).tlast);
1371         for(int i = 0; i < pDblNg->get(0); i++)
1372         {
1373             if(jroot[i])
1374             {
1375                 k++;
1376                 pDblRd->set(k, (double)i + 1);
1377             }
1378         }
1379         out.push_back(pDblRd); // rd
1380     }
1381
1382     if((meth < 3 && _iRetCount == 3) || (meth == 3 && _iRetCount == 4))
1383     {
1384         types::Double* pDblWOut = new types::Double(1,rwSize);
1385         C2F(dcopy)(&rworkSize, rwork, &one, pDblWOut->get(), &one);
1386         C2F(dcopy)(&dStructTabSize, dStructTab, &one, pDblWOut->get() + rworkSize, &one);
1387
1388         types::Double* pDblIwOut = new types::Double(1,iwSize);
1389         for(int i = 0; i < iworkSize; i++)
1390         {
1391             pDblIwOut->set(i, (double)iwork[i]);
1392         }
1393
1394         for(int i = 0; i < iStructTabSize; i++)
1395         {
1396             pDblIwOut->set(iworkSize + i, (double)iStructTab[i]);
1397         }
1398
1399         out.push_back(pDblWOut);  // w
1400         out.push_back(pDblIwOut); // iw
1401     }
1402
1403 // *** free. ***
1404     if(itol == 1 || itol == 3) // atol is scalar
1405     {
1406         free(atol);
1407     }
1408
1409     if(itol < 3) // rtol is scalar
1410     {
1411         free(rtol);
1412     }
1413
1414     free(pdYData);
1415     free(YSize);
1416     free(rwork);
1417     free(iwork);
1418
1419     if(jroot)
1420     {
1421         free(jroot);
1422     }
1423
1424     if(dStructTab)
1425     {
1426         free(dStructTab);
1427         free(iStructTab);
1428     }
1429
1430     DifferentialEquation::removeDifferentialEquationFunctions();
1431
1432     return types::Function::OK;
1433 }
1434