Multiple ways are there..
1. Fit a polynomial for the data and get the intersection points with x and y axes.
2. Use interpolation.
Eg:
clc; clear all ;
c = @(x) x.^3+x.^2-x-1 ;
N = 100 ;
x = linspace(-4,4,1000) ;
y = c(x) ;
plot(x,y) ; grid on ; hold on
% Get intersection with x-axes (y = 0)
x0 = interp1(y,x,0) ;
plot(x0,0,'*r')
% Get intersection with y-axes (x = 0)
y0 = interp1(x,y,0) ;
plot(0,y0,'*b')
In the above code, I have used inbuilt interpolation...you may write your code using for loops...
