A single Hermite cubic polynomial, fitting two points:

Clear[x1,x2,y1,y2,m1,m2,a,b]

Strategy: take a Lagrangian like approach

p[x_]:=
y1*(x-x2)^3/(x1-x2)^3+
y2*(x-x1)^3/(x2-x1)^3+
a*(x-x2)^2(x-x1)+
b*(x-x1)^2(x-x2)

Now have Mathematica solve for the a and b. What are the constraints?

p'[x1]
p'[x2]
asoln=a/.Solve[p'[x1]==m1,a]
a=asoln[[1]]
bsoln=b/.Solve[p'[x2]==m2,b]
b=bsoln[[1]]
b

Now choose a function and the abscissa to work on:

f[x_]:=Sin[x]
x1=1
x2=3
y1=N[f[x1]]
y2=N[f[x2]]
m1=N[f'[x1]]
m2=N[f'[x2]]

Let’s plot up the results:

p1=Plot[{f[x],p[x]},{x,0,4}];
p2=ListPlot[{{x1,y1},{x2,y2}},PlotMarkers->{Automatic,Medium}];
Show[p1,p2]
Plot[{f[x]-p[x]},{x,0,4}]