23 #ifndef O2SCL_ODE_BV_MULTISHOOT_H
24 #define O2SCL_ODE_BV_MULTISHOOT_H
31 #include <o2scl/astep.h>
32 #include <o2scl/astep_gsl.h>
33 #include <o2scl/ode_iv_solve.h>
34 #include <o2scl/gsl_mroot_hybrids.h>
36 #ifndef DOXYGEN_NO_O2NS
47 template <
class func_t=ode_funct<>,
48 class vec_t=ubvector,
class alloc_vec_t=ubvector,
49 class alloc_t=ubvector_alloc,
class vec_
int_t=ubvector_
int_base,
63 virtual int solve(vec_t &mesh,
int &n_func, vec_t &y_start,
64 func_t &left_b, func_t &right_b, func_t &extra_b,
65 func_t &derivs, vec_t &x_save, mat_t &y_save) {
79 int nsolve=y_start.size();
85 alloc_vec_t,alloc_t,vec_int_t> >
87 alloc_vec_t,alloc_t,vec_int_t>::
solve_fun);
91 int ret=this->
mrootp->msolve(nsolve,sx,mfm);
114 gsl_mroot_hybrids<mm_funct<> > def_mroot;
116 #ifndef DOXYGEN_INTERNAL
151 ubvector y((*this->l_n_func)),y2((*this->l_n_func));
155 for(
size_t i=0;i<(*this->
l_y_start).size();i++) {
160 for(
size_t k=0;k<(*this->
l_mesh).size()-1;k++) {
168 (*this->
l_left_b)(xa,(*this->l_n_func),sx,y);
170 for(
int i=0;i<(*this->
l_n_func);i++)
171 y[i]=sx[i+(*this->l_n_func)*(1+k)];
177 int ngrid=((*this->
l_x_save).size()-1)/((*this->l_mesh).size()-1)+1;
179 ubmatrix yysave(ngrid,(*this->l_n_func));
181 if (k!=((*this->l_mesh).size()-2)) {
182 xb=(*this->
l_mesh)[k+1]-((*this->l_mesh)[k+1]-
183 (*this->
l_mesh)[k])/ngrid;
186 this->oisp->solve_grid(xa,xb,h,(*this->l_n_func),y,ngrid,
189 for(
int i=0;i<ngrid;i++) {
190 (*this->
l_x_save)[i+k*(ngrid)]=xxsave[i];
191 for(
int j=0;j<(*this->
l_n_func);j++) {
192 (*this->
l_y_save)[i+k*(ngrid)][j]=yysave[i][j];
200 (xa,xb,h,(*this->l_n_func),y,y2,(*this->
l_derivs));
205 if (k==(*this->l_mesh).size()-2) {
206 (*this->
l_right_b)(xb,(*this->l_n_func),sx,y);
208 for(
int i=0;i<(*this->
l_n_func);i++) {
214 for(
int i=0;i<(*this->
l_n_func);i++) {
216 sy[i+k*(*this->
l_n_func)]= y2[i]-y[i];
222 (*this->
l_extra_b)(xb,(*this->l_n_func),sx,y);
223 for(
int i=0;i<(*this->
l_n_func);i++) {
224 sy[i+(int((*this->l_mesh).size()-1))*(*this->l_n_func)]=y[i];
234 #ifndef DOXYGEN_NO_O2NS