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