Numerical solution of varying coefficients ODE I have a set of observed raw data and use 2nd order ODE to fit the data 
$$y''+b_1(t)y'+b_0(t)y = 0$$
The $b_1(t)$ and $b_0(t)$ are time-dependent and I use principal differential analysis(PDA) (R-package: fda, function: pda.fd)to get the estimate of $\hat{b_1}(t)$ and $\hat{b_0}(t)$ .
To check the validity of the estimates  $\hat{b_1}(t)$ and $\hat{b_0}(t)$, I use collocation method (R-package  bvpSolve, function:bvpcol) to get the numerical solution of the ODE
$$y''+\hat{b_1}(t)y'+\hat{b_0}(t)y = 0$$
Then compare the solution from ODE above with the smoothing curve fitting of the raw data.
My question is that my numerical solution from bvpcol can caputure the shape of the fitting curve but not for the value of the function. There are different in term of some constant multiples.

See the figure of my output. The gray dot is my raw data, the red line is Fourier expansion of the raw data, the green line is numerical solution of bvpcol function and the blue line the green-line/1.62. We can see the green line can capture the shape but with values that are constant times of fourier expansion.
I fit several other data and have similar situation but different constant. I am wondering it is the problem of numerical solution of ODE or some other reasons and how to solve this problem to get a good accordance between numerical solution(green) and true Fourier expansion?
Any help and idea is appreciated!
Here is a raw data and code: 
RData is here
library(fda)
library(bvpSolve)
# load the data
load('y.RData')
tvec = 1:length(y)
tvec = (tvec-min(tvec))/(max(tvec)-min(tvec))
# create basis
fbasis = create.fourier.basis(c(0,1),nbasis=47)
bbasis = create.bspline.basis(c(0,1),norder=8,nbasis=47)
bfdPar = fdPar(bbasis)
yfd = smooth.basis(tvec,y,fbasis)$fd
yfdlist = list(yfd)
bwtlist = rep(list(bfdPar),2)
# PDA fit
bwt = pda.fd(yfdlist,bwtlist)$bwtlist
# output of estimated coefficients
beta0.fd<-bwt[[1]]$fd 
beta1.fd<-bwt[[2]]$fd 
# define the vary-coef function in terms of t
fbeta0<-function(t)eval.fd(t,beta0.fd)
fbeta1<-function(t)eval.fd(t,beta1.fd)
# define 2nd order ODE
fun2 <- function(t,y,pars) {
  with(as.list(c(y,pars)),{
    beta0 = pars[[1]];
    beta1 = pars[[2]];
    dy1 = y[2]
    dy2 = -beta1(t)*y[2]-beta0(t)*y[1]
    return(list(c(dy1,dy2)))
  })
}
# BVP 
yinit<-c(p1[1],NA)
yend<-c(p1[length(p1)],NA) 
t<-seq(tvec[1],tvec[length(tvec)],0.005)
col<-bvpcol(yini=yinit,yend=yend,x=t,func=fun2,parms=c(fbeta0,fbeta1),atol=1e-5,islin=T)
# plot output
plot(col[,1],col[,2],col='green',type='l')
points(tvec,p1,col='darkgray')
lines(yfd,col='red',lwd=2)
lines(col[,1],col[,2],col='green',type='l')
lines(col[,1],col[,2]/1.62,col='blue',type='l',lwd=2,lty=4)
legend('topleft',col=c('green','darkgray','red','blue'),
       legend=c('ODE solution','raw data','basis curve fitting','ODE solution/1.62'),lty=1)

An updated question:  are  $\hat{b_1}(t)$ and $\hat{b_0}(t)$ unique?
As I change the basis(say bspline,monomial) for expansion of $b_0(t)$ and $b_1(t)$, usually I have  very different estimated curves, but the numerical solution of ODE (green line) are very close.  That makes me believe different functional form of $\hat{b_1}(t)$ and $\hat{b_0}(t)$ can contribute the same solution. I started a new post with much more detail here.
 A: This is a linear ODE. If you have Dirichlet BC ($y(0) = y_0$, $y(1)=y_1$) the solution will be of the form $y_0*\text{sol}_0+y_1*\text{sol}_1$. How accurate are your boundary conditions compared to raw data and the basis curve fitting conditions and the boundaries?
Looking closely at the plot, you can see the ODE solution goes through the first and last points in the dataset as it should since those are used as boundary conditions. Whereas the curve fitting solution is lower on the right boundary and higher on the left boundary than the data because it is smoothing an interpolating.
If I change the bvp problem to use a value for y(1) scaled to y(1)/1.62, I get a new curve (solid black line) that overlaps the ODE solution/1.62 curve. (Note I had to take p1=y and nbasis=47 in the included code to run it.)

This is not unexpected since it is a linear ode. The real evidence is from looking at the original plot where the fit goes below the points on the right hand side. Constraining the ode solution to stop at this point is why the whole curve is 1.62 higher. If you use the values of y(0) and y(1) from the basis curve fitting, then your ODE model will match much closer. This also tells you something about the sensitivity of your model.
