Your problem is not completely defined.
There exist two ways to understand it:
 

Download A_solution.mw

@9009134

For information here is the result the statistical software R produces (fit of CP/R by a 4th degree polynomial in T)

t <- c(298., 300., 400., 500., 600., 700., 800., 900., 1000.)
cp <- c(54.187, 54.518 , 70.421 , 84.683 , 96.982 , 107.323, 116.042, 123.801 , 131.594)
cpoverR <- cp / 1.987
lm(cpoverR ~ 1+t+I(t^2)+I(t^3)+I(t^4))  # linear fit

(Intercept)            t                    I(t^2)           I(t^3)            I(t^4)  
 2.579e+00    7.263e-02    7.326e-05   -1.512e-07    6.895e-11  

Which means the model R finds is 2.579 + 7.263e-02*t + 7.326e-05*t^2 -1.512e-07*t^3 + 6.895e-11*t^4
Practically identical to the model I found "by hand".
 

A way to force LinearFit is to operate on SCALED data.
Here is an illustration for the model CP/R = f(T)
Now LinearFit and the "by hand" method (such as R) give the same result.

A remark for people from the development team that could read your thread.
I consider that a clever implementation of LinearFit should account for ill condionned problems.
Solutions are well known 

• normalize the data (as I did here)
• use methods less sensitive to bad condionning (Givens-Householder for instance)
• use othogonal polynomials and express the result in canonical basis
   

Download FITTING_after_scaling.mw

Here are fittings obtained the least squares method.

You will see that the Maple procedure LinearFit suffers a lot of problems :

1. It retuns a warning saying the (information) matrix is not of full rank.
   I show it is the case for Digits=10 (default value), but that this same matrix is of full rank for a large enough value of Digits (here 40). In my opinion this denotes a weakness of LinearFit
2. For the model S/T=f(T) (written this way, but a clever way would be to fit S=g(T) and next write S/T = g(T)/T...) LinearFit returns a completely absurd error (Error, (in Statistics:-LinearFit) the model, `S/T)`(theta), is not linear in the parameters, [])

So I completed the use of LInearFit with a "By hand" construction of the models (an extremely easy thing to do if you have basic knowledge about the ordinary least squares method).
The "rank problem" can be easily avoided and no "error" (as expected) occur when proceeding this way.

In my opinion I consider LinearFit suffers some weaknesses. 

I didn't try to use NonlinearFit instead of linearFit as a workaround because it is just stupid given all your models are linear.

Download FITTING.mw

Now, a little advice about the "good practices in linear regression".
A good way to fit statistical models is to begin by looking at the data and infer some simple model (for instance a low degree polynomial model).
If you do thys for the CP/R=f(T) model you will see it's mainly linear with only a slight discrepancy that indicates adding a quadratic term could be sufficient.
So, try to fit CP/R to the model a0+a1*T and if the result doesn't suit you fit the model a0+a1*T*a2*T^2 (and augment it if you want.

Of course (assuming no warnong or error occurs) the higher the number of parametrers (= 1+degree of the model), the lower the residual sum of squares (RSS), providing this number doesn't exceed the number of observations (9 in your case).
Let's suppose you fit polynomial models of degree d and denote by RSS(d) the value of RSS for the polynomial model of degree d ; then RSS(0) > RSS(1) > ... > RSS(8) = 0 (with your 9 data).
But be extremely careful, RSS is NOT the good criteria to use to select the best model among these d-1 models.

No other solution than the trivial one if sin(y) <> 0

PS : Of course, if y=k*Pi  (k in Z) , then any function which respects the boundary conditions is solution. for instance GenFunc(x) = (1-cos(x*2*Pi/a))
 

Download N0_Other_Solution.mw

Here is something wich looks like what you want.
Colors are randomly choosen, so the result is not always pretty (there is a randomize() command at the top of the worksheet to generate a random seed which will enable to obtain different colourings). 
You can of course define your own colors.

To obtain something closer to the picture tou present it would be necessary to rescale the radii of the vertices of each graph, a thing I didn't do here.

Left : graph1 ; middle : graph 2 ; right carteian product of g1 and g2.
Cartesian_Product_of_Graphs.mw

Many errors.

choosesystem := Maplet(.....):
Maplets[Display](maplet):

is not correct, this is probably what you wanted to do ?

maplet := Maplet(.....):
choosesystem := Maplets[Display](maplet):

Next, (look the construction of InputDialog in the help pages)

 Maplet(InputDialog['ID1']("Choose system?", 'Yes1' = Shutdown(['Yes']), 'No1' = Shutdown(['No']))):

is not correct : neither 'Ye1' nor'No1' are options of an InputDialog element. The correct syntax is (for instance, for you can exchange onapprove and oncancel)

 Maplet(InputDialog['ID1']("Choose system?", 'onapprove' = Shutdown(['Yes']), 'on cancel' = Shutdown(['No']))):

More of this : 

Shutdown(['Yes'])

means : "return the content of element named 'Yes'", which doesn't exist, so another mistake.
Mybe you wanted to return "Yes" ? In this cas write 

Shutdown(["Yes"])

Let me be more constructive now: personnaly I wouldn't have used InputDialog but RadioButton to do this.
Here is the way I would do that (but it's just my opinion and you could prefer to do it otherwise)

Download Maybe_this.mw

To complete Preben's reply,
if you only want to plot a set of points whose co-ordinates are contained in 2 vectors, possibly with undefined elements, plot or ScatterPlot can do the job without prior selection of the points with both real co-ordinates.
 

Download undefined.mw

@mlheldwein 

Here is another way using convolution and Fourier transform.
(notional examples...adjust the length of the window and the period in example 2).


 

Download Using_Convolution.mw

I don't know, maybe the expression becomes too complex with the original form of ode2 and Maple is trapped somewhere???

In such a case I always try to use a simpler form of the ode to make sure to avoid this kind of preoblem.
For instance try this

sol := dsolve({diff(Q(t),t)=A*t^B*Q(t)+C, Q(0)=S}, Q(t))

The last line replaces A, B and S by their values in the original form of ode2
 

Download workaround.mw

You will find nothing with CurveFitting ar any other smoothing strategy.
The best thing to do is to refine your "x" points (or 'i"), for instance by writting :

seqi := [seq(-10..10, h)]:  # where h << 1
for i inseqi do
...
end do:
ptsN1 := [seq(..., i in seqi)]:  # no need to use CodeTools:-Usage 

But if you do this you still won't be able to see anything for the peaks are concentrated in the ranges -4.35..-4.1 and 4.1..4.35
To verify this, just do
seqi :=  [seq(-4.35..4.1, 0.01)]:  # for the peak at negative location.

PS : Some remarks:

• There is no need to define d[i]:=i
  Replace all occurrences of d[i] by i
• There is no need to define all inner variables  under the form variable[i].
  For instance write lambda1 instead of lambda[i]
• The only exception is for f.
  Before the loop write ptsN1:=NULL:
  And inside the loop replace 
  f[d[i]]:=map(fnormal, evalf(h1[i]*(no*J1mod[i]+Omega^2*J2mod[i]))):
  by
  ptsN1 := ptsN1, [i, h1*Omega^2*J2mod ]:
  (for you have n0:=0 ; more of this map and fnormal are useless).
  After end do write 
  ptsN1 := [ptsN1]:
  This will avoid the further definition of ptsN1
• The for loop contains a lot of things that are constant and should displaced outside (a, b, ...)
  In particular 

  J2[i]:=sum(g1[i]*J1[i]+g1[i]*F2[i],k=1..1);
  Should be replaced by J2:=g1*(J1+F2); for k=..1 means k=1.

• The command 
  g1[i]:=0.5*sqrt(Pi)*tau0*exp(c2^2) *exp(-conjugate(a)*((k-1)*xr)^2)/(sqrt(conjugate(a)));

  Can be displaced outside of the loop because g1 is a constant.
  More of this, as k=1 (see above), it should be written
  g1 :=0.5*sqrt(Pi)*tau0*exp(c2^2) /(sqrt(conjugate(a)));
   
• etc
   

I do not know if this code is yours or if it is somebody's else, but I believe you must do a huge effort to (re)write it in a more clever way. 
As it is I doubt you will find someone patient enough to do this job in your place and provide you more advices than I already did.

Firstly the abscissas of your 12 points are 1, 2, ... 12, so I don't understand how you can say the period is 12, you need to add either the point x=0 or either the point x=13 to be able to say that.

Next: you say nothing about the type of function you want to draw. Dharr gave you a solution which is just one amid an infinity of solutions.
Do you wanr to interpolate the points or approximate them?

Here is a solution using interpolation.
The global function is not constructed but it is (formally) very simple to do if you need it, see here Dirac_comb
 

Download Pattern_periodisation.mw

If you prefer Dharr solution AND you 12 points, then consider the period is 11 (not 12) and modify Dharr's code this way (note n=1..5 instead of n=1..3 to obtain an interpolator, and that the period is 11. and not 11 to get a faster output).
WATCHOUT:if you do this y[1] MUST BE EQUAL to y[12]

f := (1/2)*a[0]+add(a[n]*sin(2*n*Pi*t/(11.))+b[n]*cos(2*n*Pi*t/(11.)), n = 1 .. 5);
 
fn := LeastSquares(p, t, curve = f);
fnplot := plot(fn, t = -12 .. 24, color = blue)
plots:-display(P, fnplot);

Fit_2.mw

restart:
interface(version)
Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895

with(Units[Natural]):
cos(45*Unit(degrees));
sin(90*Unit(degrees));

And in more concise way :
 

with(Units[Natural]):
d := 1*Unit(degrees): 
cos(45*d); 
sin(90*d);

Hi, 

Here is a way you maybe consider as more elegant (but I rarely use it and I prefer doing what you do yourself, that is generating uniform random numbers and "inverting" an empirical cdf...).

The UseHardwareFloats := false: instruction is usefull to convert S into a sequence of integers (look the Sample help page).
What I do not like hereare the multiple conversions this solution needs (use of round and conversion as list). I let you judge this by yourself
 

Download A_way.mw

Here is another way where I use "EmpiricalDistribution" instead of "Probability Table"
Download Another_way.mw

Try to do a few prior simplifications, for instance by setting  a = 2*h/e0, b=t/(2*(-h^2-omega^2)), and so on.
Next integrate x*expand(E) instead of x*E
Apply combine(.. trig) tio the result
...
Is it enough for you ?


 

Download Is_it_enough.mw

A solution among other 
`y'`