Effective refractive index calculation of fiber core 
*

*Can someone pls explain what is effective refractive index of fiber core and how to calculate it theoretically?

*Suppose the fiber core refractive index=1.4446 and cladding refractive index is= 1.4271, then how to  calculate the effective refractive index of the core?

*Is there any software that can directly calculate the effective refractive index of core from the core and cladding refractive index?
 A: For an intuitive answer to your first question, please see both my answer and and that of @theSkinEffect about the Effective Refractive Index. However, a more complete answer concerning the modes of a fiber is found in these lecture notes.
For your second question, what is the dimensions of the waveguide (i.e. diameter)? What is the wavelength of interest? All of these factors will determine how many and what the effective indices are for your problem. For the sake of demonstration, I will assume you are using Corning SMF-28 which has a radius of 8.3$\mu m$ and that your wavelength of interest is 1.55$\mu m$ which are very typical for the present optical communications in the telecom industry. 
To answer your question about how to calculate the effective indices, please download the free statistical software package R and drop the following code into it. I normally use Python, but since this program numerically calculates the $\beta_{pm}$, I liked the numerical solver better in R. Before you use this program, please understand better the theory of modes in a cyllindrical waveguide such as fiber because there are some assumptions made in the equations I used for the characteristics equations and you need to understand and change the order of the Bessel functions to get the correct answer.
# Written in R by Ansebbian0
# Revisions have been made to correct the mode angle
# Please be aware that the mode angle is give from the propogation direction,
# NOT from the incident angle. This can be changed by using asin instead of acos

library(pracma)


deg2rad = function(deg) {
  return((pi * deg) / 180)
}

rad2deg = function(rad) {
  return((180 * rad) / pi)
}

# Characteristic Equation for TE mode beta0
TE_eq <- function(beta0,k0,n1,n2,h){
  kapppa <- sqrt( (n2*k0)^2 - beta0^2)  # lateral wavevector
  kapppamax <- sqrt(k0**2 * (n2**2 - n1**2))
  gammma <- sqrt(kapppamax**2 - kapppa**2)
  te0 <- (besselJ(kapppa*h, 1) / (kapppa*besselJ(kapppa*h, 0))) + (besselK(gammma*h, 1, expon.scaled=FALSE) / (gammma*besselK(gammma*h, 0, expon.scaled=FALSE)))

  return(te0)
}

# Characteristic Equation for TM mode beta0
TM_eq <- function(beta0,k0,n1,n2,h){
  kapppa <- sqrt( (n2*k0)^2 - beta0^2)
  kapppamax <- sqrt(k0**2 * (n2**2 - n1**2))
  gammma <- sqrt(kapppamax**2 - kapppa**2)
  tm0 <- ((besselJ(kapppa*h, 1) / (kapppa*besselJ(kapppa*h, 0)))*(k0*n2)**2) + ((besselK(gammma*h, 1, expon.scaled=FALSE) / (gammma*besselK(gammma*h, 0, expon.scaled=FALSE)))*(k0*n1)**2)

  return(tm0)
}

# calculation of betas (set seq to small number for high resolution)
beta0_array <- function(lambda,n1,n2,h){
  k0 <<- 2*pi/lambda                        # wavenumber in free space
  beta0 <- seq(n1*k0, n2*k0, 1000)      # linspace of beta0s from smallest to largest(see pollock)
  beta0 <- beta0[-length(beta0)]            # removes the last element of beta0

  return(beta0)
}

# calculates the crossing points where there are possible zeros
pol_zeros <- function(polarization, beta0){
  intervals <- (polarization>=0)-(polarization<0)   # gives a matrix of 1 (te0>=0) and -1 (te0<0)
  differences = diff(intervals)                     # returns the difference between each variable
  izeros <- which(differences<0)                    # a neg. return from diff will mean a crossing from + to -
  X0 <- cbind(beta0[izeros],beta0[izeros+1])        # makes a Nx2 matrix of the two crossing variables

  return(X0)
}

# Numerically solves the characteristic equation [n_eff = beta0_m / k0]
modes <- function(lambda, h, n1, n2){
  beta0 <- beta0_array(lambda,n1,n2,h)
  te0 <- TE_eq(beta0,k0,n1,n2,h)                    # calls the characteristic equation for beta0 of te equation
  tm0 <- TM_eq(beta0,k0,n1,n2,h)                    # calls the characteristic equation for beta0 of tm equation

  # n_eff for TE and TM numerical calculation
  te_Xs <- pol_zeros(te0, beta0)
  nzeros_te <- dim(te_Xs)[1]

  tm_Xs <- pol_zeros(tm0, beta0)
  nzeros_tm <- dim(tm_Xs)[1]                        # finds how many rows of crossings there were

  n_eff<-matrix(data=NA, nrow=max(nzeros_te, nzeros_tm), ncol=10) # initializes a matrix for the n_eff

  colnames(n_eff) <- c('te_ne_eff', 'te_mode_beta', 'prop_dir_te_mode_angle', 'tm_ne_eff', 'tm_mode_beta', 'prop_dir_tm_mode_angle','te_kappa', 'tm_kappa', 'te_mode', 'tm_mode')

  for(i in 1:nzeros_te){                            # feeds crossing pairs into a numerical root solver
    n_eff[i,1] <- (fzero(function(x){TE_eq(x,k0,n1,n2,h)}, te_Xs[i,])$x)/k0 # fzero is a root solver in pracma
    n_eff[i,2] <- n_eff[i,1]*k0
    n_eff[i,3] <- rad2deg(acos(n_eff[i,1]/n2))
    n_eff[i,7] <- sqrt( (n2*k0)^2 - n_eff[i,2]^2)
    n_eff[i,9] <- (nzeros_te + 1) - i
  }

  for(i in 1:nzeros_tm){                        # feeds crossing pairs into a numerical root solver
    n_eff[i,4] <- (fzero(function(x){TM_eq(x,k0,n1,n2,h)}, tm_Xs[i,])$x)/k0 # fzero is a root solver in pracma
    n_eff[i,5] <- n_eff[i,4]*k0
    n_eff[i,6] <- rad2deg(acos(n_eff[i,4]/n2))
    n_eff[i,8] <- sqrt( (n2*k0)^2 - n_eff[i,5]^2)
    n_eff[i,10] <- (nzeros_tm + 1) - i
  }

  return(data.frame(n_eff)) #, nTM, TEparam, TMparam)    
}



modes(lambda=1.55e-6, h = 8.3e-6, n1=1.4271, n2=1.4446)

This is the output:

The output should be self explanatory, but if you have questions, let me know and I will tell you about how to use it. This program is not valid for the LP modes because those are composite modes. This is only valid for TE and TM. Look at Buck, Fundamentals of Optical Fibers to find the equation for the LP modes and then just replace that into the TE_eq or TM_eq function and you'll get the results.
PS: One small note. Before the above program works correctly, you need to type this in R:
install.packages("pracma")

This is the numerical solvers library which you need.
