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