Your data: data = {{0.067, 0.423}, {0.30, 0.408}, {0.60, 0.433}, {0.25, 0.3512}, {0.37, 0.4602}, {0.44, 0.413}, {0.60, 0.390}, {0.73, 0.437}, {0.8, 0.47}}; errors = {0.055, 0.0552, 0.0662, 0.0583, 0.0378, 0.080, 0.063, 0.072, 0.08}; ErrorListPlot[Transpose[{data, ErrorBar /@ errors}], PlotRange -> {0, 1}] > ![Mathematica graphics](https://i.sstatic.net/t0kVe.png) Assume that the errors are distributed Normally with standard deviation given by your errors. Then we can define diffErr[y1_, s1_, x1_, y2_, s2_, x2_] := Block[{distr, slope, error}, distr = TransformedDistribution[(v - u)/( x2 - x1), {u \[Distributed] NormalDistribution[y1, s1], v \[Distributed] NormalDistribution[y2, s2]}]; slope = Mean[distr]; error = StandardDeviation[distr]; {slope, error} ] And the slopes between points ie the numerical derivative is then slopes=Table[diffErr[data[[i, 2]], errors[[i]], data[[i, 1]], data[[i + 1, 2]], errors[[i + 1]], data[[i + 1, 1]]], {i, Length[data] - 1}] > {{-0.0643777, 0.334435}, {0.0833333, 0.287315}, {0.233714, 0.252034}, {0.908333, 0.579016}, {-0.674286, 1.26401}, {-0.14375, 0.636427}, {0.361538, 0.735933}, {0.471429, 1.53756}} ErrorListPlot[slopes, PlotRange -> {-3, 3}] > ![Mathematica graphics](https://i.sstatic.net/kY33i.png)